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1  Summary 

The  first  task  of  an  intercept  receivers  is  to  detect  the  presence  of  a  target. 
From  the  data  collected  at  mnltiple  platform  receivers  a  decision  has  to  be  made 
on  the  presence  or  absence  of  a  target.  A  common  practice  is  to  dehne  a  fnnction 
of  this  observed  data  and  compare  the  function  value  against  a  threshold  to  make 
a  decision.  This  function,  or  sometimes  also  called  as  the  test  statistic,  is  called  the 
detector.  Although  no  explicit  detector  has  been  derived  by  Fowler^,  the  approach 
taken  for  localization  suggests  that  the  detector  would  be  the  complex  ambiguity 
function  (CAF).  We  will  refer  to  this  as  the  CAF  detector  or  otherwise  simply 
as  the  time  difference  of  arrival  (TDOA)  approach  as  it  is  commonly  called.  It  is 
important  to  note  that  the  detector  that  we  are  proposing  is  also  a  TDOA  based 
detector.  What  we  are  proposing  is  a  better  technique  for  estimating  the  TDOA 
information. 

In  a  nutshell  the  CAF  detector  can  be  explained  as  follows.  From  all  the  possi¬ 
ble  receiver  combinations  only  a  small  subset  of  pairs  is  chosen.  For  the  data  from 
each  of  these  pairs  of  receivers  the  CAF  is  computed.  The  CAF  corresponding  to 
the  pair  that  yields  the  maximum  CAF  is  taken  as  the  test  statistic  and  compared 
to  a  threshold.  There  are  two  drawbacks  to  this  approach.  The  hrst  is  that  it  is 
using  only  part  of  the  data  (only  a  subset  of  all  possible  combinations).  So  the 
available  information  is  not  completely  used  which  directly  leads  to  a  poor  detec¬ 
tion  performance.  Even  if  all  the  pairs  are  used  there  is  a  second  drawback.  The 
detector  is  depending  only  on  the  correlation  factor  but  for  detection,  the  energy 
at  a  particular  receiver  is  also  an  important  factor.  This  factor  is  not  taken  into 
account  in  this  detector.  This  was  also  evident  in  some  of  our  simulations  where  a 
plain  energy  detector  outperformed  the  CAF  detector. 

Clearly,  the  currently  used  technique  is  not  processing  the  available  data  ef- 
hciently.  The  GLRT  detector  that  we  propose  processes  the  data  more  efficiently. 
We  have  shown  that  the  GLRT  detector  is  the  maximum  eigen-value  of  a  complex 
ambiguity  matrix.  This  matrix  has  the  correlation  factors  from  all  possible  com¬ 
binations  of  receiver  pairs.  Also,  the  diagonal  of  this  matrix  is  the  energy  at  each 
receiver.  So,  we  are  using  the  two  important  factors  -  energy  and  correlation,  in 
our  detector  and  so  the  performance  of  this  detector  is  better. 

Once  a  target  has  been  detected  the  next  task  of  an  interceptor  is  to  localize 
the  target.  Localization  is  the  estimation  of  the  target  position  and  velocity.  The 
TDOA  and  frequency  difference  of  arrival  (FDOA)  data  can  be  expressed  as  a  func¬ 
tion  of  the  target  position  and  velocity.  So,  if  the  TDOA  and  FDOA  information 
can  be  estimated  from  the  observed  data  then  using  these  estimates  the  target  po¬ 
sition  and  velocity  can  be  calculated.  In  the  current  technique  the  time- difference 
and  frequency-difference  values  that  maximize  the  CAF  are  taken  as  the  estimates 
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for  the  TDOA  and  FDOA  respectively.  So  for  each  of  the  chosen  pairs  of  receivers, 
corresponding  (TDOA,  FDOA)  pairs  are  estimated.  There  is  a  drawback  to  this 
approach.  For  a  given  target  location  the  (TDOA, FDOA)  pairs  corresponding  to 
each  of  the  receiver  pairs  are  all  related  and  can  take  only  a  certain  set  of  values. 
But  when  they  are  estimated  independently  this  restriction  is  not  applied  and  so 
the  estimation  is  poorer.  A  better  way  to  approach  this  is  to  express  the  OAF  as  a 
function  of  the  target  location  and  velocity.  In  our  MLE  we  avoid  these  drawbacks 
by  using  all  the  available  data  and  by  directly  estimating  the  target  location  and 
velocity  without  the  intermediate  step  of  estimating  the  TDOA  and  FDOA  values. 

It  is  often  possible  to  use  the  knowledge  of  obstacles  between  the  emitter 
and  the  sensors  and  the  azimuth  modulation  they  induce  for  the  localization  of 
the  emitter.  This  is  particularly  useful  in  urban  environments  where  there  is  a 
number  of  obstacles  such  as  buildings,  trees  etc.  These  obstacles  cause  obstruction, 
reflection,  diffraction  etc,  of  the  transmitted  signal  which  induces  some  azimuth 
modulation  in  the  signal  received  at  the  sensors.  This  kind  of  localization  is  called 
as  the  knowledge  aided  design.  We  have  derived  the  theoretical  results  that  conhrm 
the  increase  in  information  due  to  these  obstacles  under  certain  scenarios.  We  have 
also  analyzed  the  case  of  an  obstructing  obstacle. 

As  explained  previously  the  TDOA  approach  is  sub-optimal  for  localization. 
But,  the  advantage  of  the  TDOA  approach  is  that  it  requires  very  few  resources. 
For  example  since  only  a  subset  of  the  pairs  of  sensors  is  used,  data  links  are 
necessary  only  between  those  pairs.  Also,  the  TDOAs  are  estimated  at  the  local 
sensor  pairs  and  this  information  is  transmitted  to  the  central  fusion  station.  This 
requires  far  less  bandwidth  than  what  is  required  to  transmit  the  complete  signal. 
In  situations  where  the  resources  are  limited,  it  may  not  be  possible  to  implement 
the  MLE.  Under  such  circumstances,  the  TDOA  approach  is  the  only  option. 
We  have  come  up  with  an  improvement  to  the  existing  TDOA  approach  without 
requiring  any  additional  resources.  In  order  to  estimate  the  TDOA  at  a  local  pair  of 
sensors,  the  cross-correlation  function,  which  is  the  maximum  likelihood  function, 
is  maximized.  We  propose  that  the  curvature  of  the  cross-correlation  function  at 
the  peak  should  also  be  estimated  and  then  transmitted  to  the  central  fusion  station 
along  with  the  local  TDOA  estimate.  The  curvature  of  the  likelihood  function 
gives  a  measure  of  the  quality  of  the  TDOA  estimate.  Therefore,  this  curvature 
information  can  be  used  as  a  weighting  factor  when  processing  the  TDOAs  at  the 
fusion  station.  This  only  requires  one  additional  number  to  be  transmitted  to 
the  fusion  station.  We  have  found  that  the  curvature  information  is  particularly 
useful  when  some  of  the  sensors  are  operating  at  signal-to-noise  ratios  close  to  the 
break-down  range. 

In  the  following  sections  we  will  explain  each  of  the  above  mentioned  topics 
in  detail.  All  the  mathematical  derivations  are  provided  in  the  appendices. 
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2  Detection 


2.1  Introduction 

RADAR  is  an  acronym  for  radio  detection  and  ranging  with  detection  being 
the  crucial  function.  A  radar  system  illuminates  a  target  of  interest  by  transmit¬ 
ting  a  signal.  The  echoes  that  are  reflected  back  are  fed  to  a  detector  which  makes 
a  decision  on  the  presence  or  absence  of  a  target.  These  radars  are  called  active 
radars.  There  is  another  class  of  radars  called  passive  radars  or  interceptors  which 
silently  listen  for  transmissions  from  other  active  devices  such  as  active  radars, 
jammers,  beacons,  etc.  Intercept  receivers  are  most  desirable  in  hostile  situations 
due  to  their  covert  nature  [1].  In  addition  there  is  only  a  one-way  power  loss  at 
the  intercept  receiver  as  compared  to  the  two-way  loss  at  the  transmitting  active 
radar.  On  the  other  hand,  since  the  angle  of  signal  arrival  and  the  signal  itself 
are  generally  unknown  to  the  interceptor,  efficient  processing  techniques  such  as 
matched  hltering  cannot  be  implemented.  Also,  many  modern  active  radar  systems 
are  designed  with  low  probability  of  intercept  (LPI)  features.  They  incorporate 
physical  attributes  such  as  frequency  variability,  infrequent  scanning,  etc.  to  re¬ 
duce  the  probability  of  interception  by  an  interceptor  and  signal  design  attributes 
such  as  low  power,  wide  bandwidth,  etc.  to  decrease  the  probability  of  detection 
and  parameter  identification  at  the  interceptor  [2].  This  leads  to  the  need  for 
implementing  highly  efficient  detectors  in  the  interceptor  systems. 

A  simple  radar  system  with  a  single  transmitter  and  a  single  receiver  both  at 
the  same  physical  location  is  called  a  monostatic  radar.  In  general,  its  performance 
is  inferior  to  a  multistatic  radar  system.  An  active  multistatic  radar  system  has  one 
or  more  transmitters  and  many  spatially  separated  receiving  stations.  A  passive 
multistatic  radar  consists  of  only  a  network  of  distributed  sensors.  Such  a  system 
of  multiple  receive  platforms  has  a  higher  probability  of  intercepting  the  signals  of 
interest.  Also,  multiple  platforms  provide  more  data  samples  over  a  given  interval 
of  time  which  increases  the  probability  of  detection.  Each  of  these  receive  plat¬ 
forms  can  perform  some  kind  of  processing  on  the  received  signal.  This  is  called 
decentralized  target  detection  [3,  4].  These  processed  signals  from  the  individual 
receive  platforms  are  communicated  to  a  central  processor  where  they  are  further 
processed  to  arrive  at  a  global  decision.  Quite  often  the  local  receive  platforms 
process  the  received  signal  to  arrive  at  a  local  decision.  These  local  decisions  are 
communicated  to  a  central  processor  for  decision  fusion.  Several  such  decision 
fusion  techniques  have  been  proposed  and  used  over  the  years  [5,  6,  7].  Despite 
some  practical  advantages  such  as  requiring  low  bandwidth  data  links  between 
receive  platforms  and  less  processing  power,  the  decentralized  detection  approach 
has  an  obvious  performance  loss  due  to  the  absence  of  cross-platform  correlation 
information.  On  the  other  hand,  centralized  detection  has  better  performance 
but  requires  more  resources.  Large  bandwidth  data  links  are  required  to  trans¬ 
fer  the  received  signals  from  all  the  receiving  stations  to  the  fusion  center.  High 
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speed  signal  processors  are  required  to  process  the  large  amount  of  data  available 
at  the  fusion  center  for  real-time  detection.  Since  centralized  detection  relies  on 
the  cross-platform  correlation  information,  the  receive  platforms  must  be  precisely 
synchronized  in  time.  In  [8],  a  centralized  detector  for  a  particular  model  is  given. 
Here  we  use  a  model  that  is  similar  to  the  model  used  by  Stein  in  [9]. 

When  a  signal  is  transmitted  by  an  active  device,  from  here  on  referred  to 
as  the  target,  it  reaches  each  of  the  receiving  stations  after  a  certain  amount  of 
time,  which  is  the  propagation  delay.  If  the  target  is  moving  with  respect  to  the 
receiving  station,  then  there  is  a  Doppler  shift  in  the  frequency  of  the  signal.  Also, 
the  amplitude  and  phase  of  the  signal  are  changed  due  to  propagation  through  the 
channel.  Taking  all  these  factors  into  account  and  assuming  there  is  no  noise,  the 
complex  envelope  of  the  signal  that  is  received  at  the  ith  receiving  station,  after 
sampling,  is 


Ais[n  —  rij]  exp 


/ j2nki{n  -  Ui) 

V  N 


where  s[n]  is  the  signal  emitted  by  the  target,  Ai  =  Aie^"^*  is  the  change  in  am¬ 
plitude  and  phase,  rii  is  the  propagation  delay  and  ki  is  the  Doppler  shift.  (For 
simplicity  we  assume  a  discrete-frequency  Doppler  shift  of  for  some  large  N .) 
Note  that  we  are  using  to  represent  complex  variables.  An  important  differ¬ 
ence  between  detectors  in  active  and  passive  radars  is  that  in  the  active  case  the 
received  signal  is  an  echo  of  the  transmitted  signal.  Thus,  the  signal  can  be  mod¬ 
eled  as  known  with  unknown  parameters.  In  the  passive  detection  case,  however, 
the  signal  itself  is  unknown.  In  some  cases  the  signal  is  modeled  as  a  stochastic 
process  with  unknown  parameters  [10].  We  model  the  signal  as  deterministic  and 
completely  unknown,  i.e,  our  signal  modeling  assumptions  coincide  with  those  of 
Stein  [9]. 

In  the  next  section  we  give  a  detailed  description  of  the  problem  and  model  it 
as  a  statistical  hypothesis  test.  In  Section  2.3  we  hnd  the  GLRT  for  the  described 
hypothesis  test.  We  compute  the  GLRT  analytically  for  the  special  case  of  2  sen¬ 
sors.  Here  we  show  that  the  maximum  likelihood  estimate  (MLE)  for  the  delay  and 
Doppler  is  obtained  by  maximizing  the  GAF.  In  [9],  Stein  addressed  the  problem 
of  differential  delay  and  Doppler  estimation  for  the  case  of  two  sensors.  He  also 
arrived  at  the  same  result  that  the  MLE  for  the  differential  delay  and  Doppler  is 
obtained  by  maximizing  the  GAF.  He  did  not  address  the  problem  for  more  than 
two  sensors.  Also,  he  did  not  address  the  problem  of  target  detection.  In  Section 
2.4  we  compare  the  performance  of  the  GLRT  against  some  commonly  used  detec¬ 
tors.  In  Section  2.5  we  make  further  assumptions  to  simplify  the  problem  model 
and  derive  the  respective  GLRTs.  We  provide  conclusions  in  Section  2.6. 


2.2  Methods,  Assumptions,  and  Procedures 

The  problem  we  are  addressing  can  be  described  as  follows.  We  have  M 
intercepting  sensors  placed  at  multiple  locations.  Each  of  these  sensors  collects  N 
time  samples  in  a  given  interval  of  time.  The  total  MN  samples  that  are  collected 
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are  available  for  processing  at  a  fusion  center.  As  is  usual  in  practice,  we  assume 
N  >  M.  We  further  assume  that  the  noise  at  each  of  the  receiver  stations  is 
white  Gaussian  and  also  that  the  noise  at  a  receiver  station  is  independent  from 
noise  at  all  the  other  receiver  stations.  For  simplicity  we  assume  that  the  variance 
of  noise  is  =  1.  In  situations  where  the  noise  does  not  meet  these  criteria, 
the  observations  can  be  pre-whitened  [11]  before  processing  so  that  the  above 
assumptions  hold.  Now,  the  task  of  the  detector  at  the  fusion  center  is  to  detect 
the  presence  of  an  unknown  signal  in  the  MN  observed  samples. 

Let  s[n],n  =  0, 1  ■  •  •  ,  A^  — 1  be  the  nth  complex  time  sample  of  the  transmitted 
signal.  A  time  delayed  and  Doppler  shifted  version  of  this  signal  and  with  a 
different  amplitude  and  phase  reaches  the  sensors.  Let  be  the  change 

in  the  amplitude  and  phase.  Let  n*  be  the  discrete  time  delay  and  ki  be  the  discrete 
Doppler  shift  .  Let  Wi[n]  be  the  nth  complex  time  sample  of  the  additive  noise  at 
the  Ah  receiver  station.  If  rj[n]  is  the  nth  time  sample  of  the  observation  at  the 
Ah  receiver  station  then  we  can  write 

ri[n]  =  Ais[n  -  n^j  exp  j  yj.[n]  n  =  0, 1,  •  •  •  ,  -  1 

i  =  0,1,---  ,M-  1. 

When  there  is  no  signal,  the  observation  is  just  noise,  rj[n]  =  Wi[n\.  So,  the 
hypothesis  test  for  the  detection  problem  can  be  written  as 

Ho  :  fi[n]  =  Wi[n] 

Hi  :  rjn]  =  Ais[n  —  n*]  exp  j  _|_  yjjn]  n  =  0, 1,  •  •  •  ,  A^  —  1 

i  =  0,l,---  ,M-1. 

Notice  that  here  we  are  modeling  the  received  signal  as  unknown  but  deter¬ 
ministic  with  unknown  parameters  Ai,ni,ki,i  =  0,1,--  -  ,M  —  1  and  s[n],n  = 
0, 1,  •  •  •  ,  A^  —  1.  Now,  if  we  let  u  =  exp(^)  and  W  be  the  N  x  N  matrix  W  = 
diag(a;°,  and  let  P  be  an  A^  x  A^  permutation  matrix  dehned  as 

\P]ij  =  lii  i  =  j  +  1  and  0  otherwise,  i  =  0, 1,  •  •  •  ,  A^  —  1,  j  =  0, 1,  •  •  •  ,  A^  —  1  and 
[P]o,7v-i  =  1,  then  we  can  write  the  hypothesis  test  in  vector  form  as  follows 

■Ho ;  s  =  0 
Hi  ;  s  7^  0 


where 

ii  =  AiP^^w'^^s  +  Wi  i  =  0,1,---  ,M- 1 

and  fi  =  [  ri[0]  A[l]  •••  ri[A^  -  1]  s  =  [  s[0]  s[l]  ■■■  s[A^  -  1]  and 

Wj  =  [  tCi[0]  hi^l]  "D^A^  —  1]  ]  •  The  permutation  matrix  P  circularly 

shifts  s.  This  causes  an  effect  as  if  s  was  periodic  with  period  equal  to  N  samples. 
So,  as  in  [9],  we  assume  that  the  signal  s  is  non-zero  only  in  an  interval  that  is 
smaller  than  N  samples  and  that  the  discrete-time  delays  are  relatively  small  com¬ 
pared  to  N.  Now  let  f  =  [  f Q  w  =  [  Wg  wf  ■  ■  ■  ]^, 
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A  =  [  Ao  Ai  •••  Am-1  n  =  [  no  ni  •••  n^-i 

k  =  [  /co  ki  ■  ■  ■  kM-i  ]'^  and 


T 


H(A,  n,  k) 


The  hypothesis  test  can  be  written  as  follows 


T-Lo  :  s  =  0 
Hi  ;  s  ^  0 


where 


r  =  H(A,  n,  k)s  +  w 


and  w  has  a  complex  normal  distribution  with  zero  mean  and  the  identity  matrix 
as  the  covariance  matrix,  i.e,  w  ~  CA/’(0,Imv),  ^mn  is  an  MN  x  MN  identity 
matrix.  Here  A  is  M  x  1,  s  is  A  x  1,  n  is  M  x  1,  k  is  M  x  1,  and  are  all  assumed 
unknown. 


2.3  GLRT  Detector 

In  hypothesis  testing  problems  where  the  probability  density  functions  (PDFs) 
under  both  the  hypotheses  are  completely  known  the  Neyman-Pearson  (NP)  detec¬ 
tor  is  the  uniformly  most  powerful  (UMP)  detector  [12].  When  there  are  unknown 
parameters  in  the  PDFs,  the  performance  of  the  NP  detector  depends  on  the  true 
value  of  these  parameters.  There  are  two  common  ways  to  deal  with  these  un¬ 
known  parameters.  They  can  be  modeled  as  random  variables  with  some  PDF 
and  then  integrated  out  or  they  can  be  modeled  as  unknown  but  deterministic  and 
replaced  with  their  MLE.  In  the  GLRT  the  unknown  parameters  are  replaced  by 
their  MLEs  under  the  different  hypotheses.  Asymptotically,  the  GLRT  is  the  UMP 
test  among  all  tests  that  are  invariant  [12] .  The  likelihood  ratio  for  the  previously 
described  hypotheses  test  is  given  by 


L(f) 


p(r;s,  A,n,k,?^i) 


(1) 


If  we  replace  the  unknown  parameters  A,  s,  n  and  k  with  their  respective  MLEs 
A,  s,  h  and  k  then  the  GLRT  decides  Hi  if 


p(f;s,  A,h,k,?^i) 

As  derived  in  Appendix  A,  the  GLRT  test  statistic  is 

In  Loir)  =  maxAmax(B(n,  k)) 

n.k 


(2) 

(3) 
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where  Amax  is  the  maximum  eigenvalue  and  B  is  the  M  x  M  complex  cross¬ 
ambiguity  matrix  (CAM)  given  by 


[B]ij  = 


(4) 


where  H  is  conjugate  transpose.  Since  B  is  positive  dehnite,  the  maximum  eigen¬ 
value  is  real  and  positive.  Therefore  the  GLRT  decides  that  a  signal  is  present  if 
the  maximum  eigenvalue  of  the  CAM,  when  also  maximized  over  time  delay  and 
Doppler,  is  greater  than  a  threshold.  It  should  be  noticed  that  the  diagonal  ele¬ 
ments  in  the  CAM  are  the  energy  terms  at  each  of  the  receivers,  i.e,  and  the 
off-diagonal  terms  are  the  CAF  values  between  all  the  pairs  of  sensors.  For  two 
hnite  length  discrete  time  complex  signals  ro[nj^and  ri[n],  n  =  0, 1,  •  •  •  ,N  —  1  the 
energy  in  each  of  the  signals  is  given  by  Si  =  J2n=o  IdN  P)  *  =  0, 1  respectively  and 
the  CAF  is  a  two  dimensional  function  of  differential  delay  (An)  and  differential 
Doppler  shift  (A/c)  between  the  signals  and  is  dehned  as 


CAF{An,Ak)  =  ^  ro[n]rj['[n -I- An]  exp 

where  *  represents  complex  conjugate.  It  is  important  to  note  that  the  CAM 
contains  CAF  terms  for  all  possible  sensor  pair  combinations  and  not  just  a  selected 
set  of  pairs.  So,  the  CAM  can  also  be  written  as  [B]jj  =  Si  ii  i  =  j  and  [B]jj  = 
CAFij  if  i  ^  j  where  CAFij  is  the  CAF  of  the  observations  at  sensor  i  and  sensor 
j.  Another  interesting  result  is  that  when  the  correlation  information  is  zero,  all 
the  off  diagonal  terms  become  zero  and  so  the  maximum  eigenvalue  is  simply  the 
maximum  of  the  energies  of  the  sensors.  Hence  the  GLRT  simply  reduces  to  a  type 
of  energy  detector. 

At  hrst  glance,  the  maximization  of  the  eigenvalue  appears  to  be  on  a  2M 
dimensional  space.  This  computation  can  be  prohibitive  when  the  number  of 
sensors  is  large.  But  it  should  be  noted  that  the  true  values  of  the  delay  and 
Doppler  he  in  a  much  smaller  space.  This  is  explained  by  the  relation  of  the 
delay  and  Doppler  to  the  target  location  and  velocity.  If  we  assume  that  the 
location  and  velocity  of  the  sensors  is  known,  then  the  delays  to  the  sensors  are  a 
function  of  the  target  location  and  the  Doppler  shifts  are  a  function  of  the  target 
location  and  velocity.  So,  if  (xt,  Vt,  zt)  are  the  three  dimensional  coordinates  of 
the  target  location  and  {vx,Vy,Vz)  are  the  target  velocity  components  in  the  x,y 
and  2;  directions  respectively,  then  we  can  write  the  delays  as  n{xT,  Vt:  zt)  and  the 
Doppler  shifts  as  k{xT,yT,  ZT,Vx,Vy,Vz).  Putting  these  back  in  equation  (3)  we 
have 

InLG(f)  =  max  \^s,^(B{xT,yT,  ZT,Vx,Vy,Vz))  (6) 

XT,yT,ZT,Vx,Vy,Vz 

Hence,  the  maximization  is  at  most  on  a  six  dimensional  space  and  can  be  per¬ 
formed  numerically.  Any  additional  information  about  target  location  and  velocity, 
viz.  target  is  on  the  ground  or  target  is  stationary,  further  reduces  this  space.  No¬ 
tice  that  the  (f^,  j/r,  zt,  Vx,  Vy,  Vz)  that  maximize  the  test  statistic  are  the  MLEs  of 
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the  target  position  and  velocity.  Therefore,  we  have  simultaneously  estimated  the 
target  location  and  velocity  as  well.  We  will  investigate  the  localization  problem 
in  section  3.  Weiss  investigated  the  problem  of  localization  of  narrowband  radio 
frequency  transmitters  in  [13]  and  arrived  at  a  similar  result. 

The  fact  that  the  GLRT  test  statistic  is  a  function  of  the  maximum  eigenvalue 
of  the  CAM  can  be  further  understood  as  follows.  The  problem  can  be  viewed  as  a 
rank  one  approximation  of  an  observation  matrix.  Let  R  be  the  N  xM  observation 
matrix  given  by  R  =  [  fg  fi  ■  ■  •  tm-i  ]•  Assuming  there  are  zero  delays  and 
Doppler  shifts,  we  have 

R  =  sA^  +  Wtv  (7) 

where  W^r  is  the  N  x  M  noise  matrix  given  by  W^r  =  [  wq  wi  ■  ■  ■  wm-i  ]• 

_  ~  T 

Note  that  the  signal  component  sA  is  a  rank  one  matrix  with  one  singular  value 

=  \l  s^sA^A  =  so  that  the  maximum  eigenvalue  is  Amax  =  s^sA  A. 

The  noise  W at  causes  R  to  be  full  rank.  Hence,  the  detector  attempts  to  extract 
a  characteristic  of  the  signal.  It  is  conjectured  that  in  the  presence  of  multiple 
targets,  say  n  targets,  the  test  statistics  will  be  the  maximum  n  eigenvalues  of  the 
CAM. 


2.3.1  Example:  A  simple  2-sensor  case 

We  have  seen  that  the  CLRT  statistic  is  a  function  of  the  maximum  eigenvalue. 
In  general,  it  is  difficult  to  compute  the  eigenvalues  analytically,  but  for  the  special 
case  of  2-sensors  it  is  possible.  When  M  =  2  the  CAM  is  a  2  x  2  matrix  and  so 
the  eigenvalues  are  computed  as  follows.  We  have 

1  r  -| 

Since  (note  that  P^P  =  Iat)  and  =  Iat,  where  Iat  is 

an  A  X  A  identity  matrix,  we  have 

^  ^  [  f"f„  r5'P”«w‘"(w'=‘)"(P”-)»fi  ■ 

’  [  ffP"‘w''‘(w'*)''(P"«)"fo  rffi 

Notice  that  the  diagonal  elements  are  the  energies  at  individual  sensors  and  the 
off  diagonal  terms  are  the  cross  ambiguity  terms.  Since  B(n,  k)  is  a  2  x  2  matrix 
it  has  two  eigenvalues  which  can  be  determined  analytically.  The  maximum  of  the 
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two  real  and  positive  eigenvalues  is 


f^fo  +  ff  fi  +  y [f^fo  -  +  4 


Amax(B(n,  k))  = 

In  order  to  maximize  Ajnax(B(n,  k))  over  (n,  k),  we  have  to  maximize 


(8) 


If  we  let  Ak  =  ki  —  ko  and  An  =  ni  —  Uq  then,  since 

we  have 


max 

n,k 


rf  P^i  (P«o)/^fo 


max 

An,Ak 


pK-nP 


N-1 

fo[n]fj[‘[n  +  An]  exp 

n=0 

(9) 

The  above  expression  is  the  magnitude  of  the  CAF  in  discrete  time  and  the  An 
and  Ak  that  maximize  it  are  the  MLEs  of  An  and  Ak  respectively.  A  similar 
result  was  derived  in  [9]  for  the  2  sensor  case  in  continuous  time,  but  the  CAF  was 
maximized  only  to  obtain  the  MLEs  of  An  and  Ak  and  the  detection  problem  was 
not  addressed.  In  [14]  however.  Holt  derived  a  similar  expression  for  the  GLRT.  He 
showed  that  when  the  signal  is  modeled  as  deterministic  and  completely  unknown, 
the  GLRT  test  statistic  is  a  weighted  sum  of  the  energies  in  the  observations  at  each 
sensor  plus  the  real  part  of  the  appropriately  weighted  cross  ambiguity  function. 
Here  we  have  derived  the  GLRT  test  statistic  to  be,  from  (8)  and(9) 


/  j2'KAkn\ 
— 


=  max 
An,Ak 


In  Loir) 


r^rp  +  rf  ri 
2 


^0^0  -  rf  i-i 
2 


2 


+  max  \CAF{An,Ak)f 

An,Ak 


(10) 


where  CAF{An,  Ak)  =  J2n=o  +  An]ro[n]exp  a^d  is  the  CAF  pre¬ 

viously  given  in  (5).  Therefore,  the  GLRT  incorporates  the  CAF  as  well  as  the 
energy  information. 


2.4  Results  and  Discussion 

We  compared  the  performance  of  the  GLRT  against  two  commonly  used  de¬ 
tectors.  One  detector  computes  the  energy  of  the  observations  individually  at  each 
sensor  and  when  the  energy  at  any  one  sensor  exceeds  a  predetermined  threshold, 
the  target  is  declared  as  present.  In  mathematical  terms,  if  £*,  i  =  0, 1,  •  •  •  ,  M  —  1 
are  the  energies  of  the  observations  at  the  M  sensors,  then  the  detector  decides 
that  a  target  is  present  if  max{To,  Si,  -  ■  ■  ,  Sm-i}  >  Is-  This  is  called  the  maximum 
energy  detector.  It  can  be  noticed  that  this  detector  is  based  solely  on  the  energy 
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information.  The  other  detector  is  one  that  is  based  solely  on  the  cross-sensor  cor¬ 
relation  information.  Here  a  sensor  is  hxed  as  the  reference  sensor  and  the  CAFs 
are  computed  between  observations  from  the  reference  sensor  and  all  the  other 
sensors.  When  the  maximum  magnitude  over  possible  delay  and  Doppler  shifts  of 
any  one  of  these  CAFs  exceeds  a  threshold,  the  target  is  declared  as  present.  In 
mathematical  terms,  assume  sensor  0  is  the  reference  sensor  and  An*  =  (uj  —  uq) 
and  A/cj  =  {ki  —  fco)  are  the  difference  in  delay  and  Doppler  at  sensor  i  and  sensor 
0  respectively.  Then  if 


|CAFi|  =  max 

Arii,  Afci 


N-1 


^  ro[n]r*  [n  +  Arij]  exp 


n=0 


^j27rAfcjn 


N 


,  z  =  l,2,--  -  ,M-1 


are  the  maximum  magnitudes  over  delay  and  Doppler  of  the  CAFs  between  sensor 
0  and  all  the  other  sensors,  then  a  target  is  declared  as  present  if 


max{|CAFi|,  ICAF2I,  •  •  •  ,  \CAFm-i\}  >  Icaf- 
This  is  referred  to  as  pair-wise  maximum  CAF  detector. 


Physical  setup  of  the  sensors  and  the  target 


Figure  1.  Physical  placement  of  the  sensors  and  the  target  position  used  for  sim¬ 
ulation. 

For  the  purpose  of  simulation  a  set  of  11  sensors  were  placed  in  a  configuration 
as  shown  in  Figure  27.  To  simplify  the  computations  we  assumed  that  the  target 
is  stationary  and  so  the  Doppler  shift  is  zero.  However,  delays  are  incorporated. 
We  used  a  Gaussian  pulse  for  the  signal  which  is  shown  in  Figure  12.  The  length 
of  the  signal  in  time  is  10  /is  and  its  bandwidth  is  approximately  0.5  MHz.  The 
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Unknown  Transmitted  signal 


Figure  2.  A  Gaussian  pulse  that  is  used  as  the  unknown  transmitted  signal 

signal  was  sampled  at  a  rate  of  =  3  MHz  to  collect  30  non-zero  samples.  A 
total  observation  interval  of  A^  =  450  samples  was  required  in  order  to  allow  the 
signal  to  reach  all  the  sensors.  This  is  shown  in  Figure  9.  White  Gaussian  noise 


Signals  reaching  the  sensors  -  no  noise  case 


Figure  3.  Signal  reaching  different  sensors  at  different  times  with  different  atten¬ 
uations. 

was  used  as  the  additive  noise  at  the  sensors.  The  noise  variance  was  adjusted  so 
that  the  average  engery-to-noise  ratio  (AENR)  is  10  dB.  The  AENR  is  the  ratio 
of  the  energy  in  the  signal  to  the  noise  power  at  each  sensor,  averaged  over  all 
the  sensors,  i.e.,  if  Esi  =  energy  of  the  signal  at  the  ith 

sensor  and  is  the  noise  variance  at  the  Ah  sensor,  then  the  signal  to  noise  ratio 

averaged  over  M  sensors  is  given  by  10 log  J2iLo^  fs-j-  •  Here  we  assumed  a 

single  target  is  located  at  (130,75)  km  and  the  maximization  of  the  GLRT  statistic 
was  done  over  {x,y).  We  used  a  grid  search  to  maximize  the  test  statistic.  The 
comparison  receiver  operating  characteristics  (ROG)  curves  are  shown  in  Figure 
4.  We  ran  1000  simulations  for  each  of  the  detectors  to  generate  the  ROG  curves. 
The  grid  search  was  performed  over  a  grid  of  size  3  km  x  3  km  around  the  true 
target  location  with  a  grid-point  distance  of  62.5  m.  The  computation  times  in 
MATLAB  for  the  pair-wise  maximum  GAF  detector,  maximum  energy  detector 
and  the  GLRT  are  given  in  Table  2.4. 
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1 


ROC  curves  for  average  ENR  =  10dB  per  sensor 


Figure  4.  Comparison  of  GLRT  against  the  maximum  energy  detector  and  the 
pair-wise  maximum  CAF  detector 


Table  1.  Comparison  of  the  computation  times  for  the  three  detectors. 


pair-wise  maximum  GAF  detector 

maximum  energy  detector 

GLRT 

7.18  sec 

1.58  sec 

778.78  sec 

Incorporation  of  the  Doppler  will  increase  the  computation  time  further  and  hence 
has  not  been  included  at  this  time.  Usage  of  efficient  maximization  techniques 
may  reduce  the  computation  time  at  which  point  the  Doppler  parameter  could 
be  included.  At  this  point  we  have  investigated  the  GLRT  itself  and  have  left 
the  details  of  its  efficient  numerical  computation  for  a  future  paper.  It  can  be 
noticed  that  the  performance  of  the  GLRT  is  very  much  better  than  either  of 
the  two  commonly  used  detectors.  For  a  probability  of  false  alarm  PpA  =  0.01 
the  probability  of  detection  Pq  is  about  0.05  for  the  maximum  energy  detector 
and  0.02  for  the  pair-wise  maximum  GAF  detector  but  is  0.45  for  the  GLRT. 
Therefore,  the  GLRT  which  is  a  function  of  the  energy  at  each  sensor  combined 
with  the  cross-sensor  correlation  is  better  than  detectors  that  are  solely  based  on 
either  energy  or  cross-sensor  correlation.  Altes  arrives  at  the  same  conclusion  in 
[15].  His  results  indicate  that  cross-correlation  alone  is  usually  sub  optimum.  He 
concludes  that  detection  should  use  a  weighted  sum  of  pairwise  cross-correlation 
between  subarrays  and  energy  detection  at  each  subarray. 

Next,  in  Figure  5,  the  same  emitter  signal  was  used  for  the  case  of  M  =  2 
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ROC  curves  for  average  ENR  =  lOdB  per  sensor 


Figure  5.  GLRT  vs  maximum  energy  detector  for  the  two  cases  of  M=2  and  M=ll. 

and  compared  against  the  case  of  M  =  11.  It  should  be  noticed  that  when  M  =  2 
the  GLRT  does  only  slightly  better  than  the  energy  detector  as  there  is  very  little 
cross-correlation  information  but  when  M  is  increased  to  11  the  performance  of 
GLRT  exceeds  that  of  energy  detector  considerably.  The  AENR  at  each  sensor  was 
set  to  10  dB.  In  real  world,  as  the  signal  travels  from  the  target  to  the  sensor  the 
signal  power  is  attenuated,  which  is  called  the  propagation  loss.  In  our  simulation, 
while  generating  the  observations  at  each  sensor,  we  accounted  for  this  propagation 
loss  also  when  computing  the  AjS  by  making  the  RjS  inversely  proportional  to  the 
distance. 

2.5  Some  Simpler  Models 

The  GLRT  detector  that  was  previously  derived  is  a  complete  solution  ac¬ 
counting  for  time  delays  and  Doppler  shifts.  Depending  on  the  specihc  problem 
appropriate  assumptions  can  be  made  in  the  GLRT  to  arrive  at  simpler  models 
that  can  be  more  easily  implemented. 

2.5.1  Simple  Bilinear  Model 

Assume  a  situation  where  the  target  and  the  sensors  are  stationary  -  hence 
there  is  no  Doppler.  Also,  assume  that  if  a  target  is  present,  its  location  is  a  priori 
known.  In  such  cases  the  time  delay  and  Doppler  parameters  can  be  dropped  and 
a  simple  bilinear  (multiplicative)  model  [16]  can  be  used.  For  such  a  model  the 
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signal  model  and  hypothesis  test  are  given  by 

f  =  H(A)s  +  w 
■Ho ;  s  =  0 
Hi  ;  s  7^  0 

where  H(A)  =  [  AqIat  Ailjv  •  •  ■  Am-i^-n  ]•  Dropping  the  delay  and  Doppler 
terms  in  equation  (3),  we  have  the  GLRT  test  statistic  for  this  model  as 

InLG(f)  =  Amax(B) 

where  Amax  is  the  maximum  eigenvalue  and  B  is  cross-sensor  correlation  matrix 
given  by 

*'o  ^0  ^0  '  '  '  ^0  1 

ff  fo  ff  fi  •  ■  ■  ff  fM-1 

^m-Ai  ■  ■  ■  ^m-Am-1 

Notice  that  the  principal  diagonal  elements  of  B  are  the  energies  at  each  of  the 
sensors.  The  off-diagonal  elements  are  the  cross-correlation  terms.  These  results 
are  analogous  to  the  results  in  Section  2.3  except  that  here  B  is  not  a  function  of 
the  delay  and  Doppler. 


2.5.2  Classical  Linear  Model 

This  assumes  furthermore  that  the  signal  arriving  at  the  sensors  has  the  same 
amplitude  and  phase  which  we  incorporate  into  the  unknown  signal  s.  So,  assuming 
Ai  =  1  for  i  =  -  ,M  —  1,  we  have  the  classical  linear  model  as 

f  =  Hs  -I-  w 


where  the  MN  x  N  matrix  H  =  [  I^r  Iv  ■  ■  ■  Iv  is  known  and  f,  s,  w  are 
same  as  dehned  in  Section  2.2.  The  hypothesis  test  for  this  model  is  given  by 

Ho  :s  =  0 
Hi  ;  s  7^  0 

The  GLRT  for  this  hypothesis  is  to  decide  Hi  if  [11] 


T(f) 


172 


where 


^  M-l 
i=0 


(11) 
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is  the  MLE  of  s  under  "Hi  (  note  that  =  MIjv).  Combining  the  two  equations, 

we  have 


T(f)  = 


M-l 

^i«f. 

i=0 


ri[n\s  [n\ 


(12) 


1/2  1/2 

which  is  an  estimator-correlator  summed  over  all  the  sensors.  When  the  additive 
noise  is  Gaussian,  the  T(f)  has  a  central  chi-squared  distribution  under  "Hq  and  a 
noncentral  chi-squared  distribution  under  "Hi.  The  exact  detection  performance  is 
given  by 

Pfa  = 


(Y) 


(13) 


where  <5^2^  is  the  right-tail  probability  of  a  random  variable  with  central  chi- 
squared  distribution  with  2N  degrees  of  freedom  and  Q^,2  is  the  right-tail  prob¬ 
ability  of  a  random  variable  with  noncentral  chi-squared  distribution  with  2N 
degrees  of  freedom  and  with  A  as  the  noncentrality  parameter.  The  noncentrality 
parameter  is 

A  =  s^(H^H)s  =  Ms^s  (14) 

which  is  the  total  signal  energy  at  all  the  sensors. 


2.5.3  Total  Energy  Detector 

In  Section  2.4  we  have  seen  the  maximum  energy  detector  that  is  commonly 
used  where  the  energies  at  each  of  the  sensors  are  independently  compared  against 
respective  thresholds.  This  is  decentralized  detection.  If  the  energies  computed 
at  each  sensor  are  all  summed  at  a  fusion  center  and  used  as  a  test  statistic,  then 
we  have  the  total  energy  detector.  The  fact  that  the  signal  received  at  each  of 
the  sensors  is  originating  from  the  same  source  is  ignored  here  also  in  the  detector 
design.  If  Sj  =  [  Si[0]  Si[l]  •  •  •  —  1]  ]  is  the  source  signal  at  the  ith.  sensor 

[rri  rrf  rri  -i  'J^  ' 

Sq  s(  ■  ■  ■  J  ,  then  using  the  same  dehnitions  for  r  and  w  as 

in  Section  2.2,  we  have  a  classical  linear  model  for  this  problem  as 

f  =  s  -|-  w 


Notice  that  here,  unlike  the  previous  case,  we  have  modeled  the  source  signal  for 
each  sensor  as  a  different  unknown  parameter.  The  hypothesis  test  for  this  model 
can  be  written  as 

Pq  ;  s  =  0 
Pi-.s^Q 

The  GLRT  for  this  hypothesis  test  is  to  decide  "Hi  if  [11] 


T(f) 


5^5 

172 


>7 
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where  s 
have 


r  is  the  MLE  of  s  under  "Hi.  Putting  this  in  the  above  equation,  we 


T(f) 


1/2 


M-1  N-1 


1/2 


(15) 


and  should  be  compared  to  (12).  When  the  additive  noise  is  Gaussian,  here  also 
the  T(f)  has  a  central  chi-squared  distribution  under  Hq  and  a  noncentral  chi- 
squared  distribution  under  "Hi.  But  the  degrees  of  freedom  here  is  2MN.  The 
exact  detection  performance  of  this  detector  is  given  by 


PpA  — 

Pd  =  Q 


^X2MN 


(7' 


X2MnW 


(7 


(16) 


where  the  noncentrality  parameter  is 


A 


1/2 


A  comparison  of  this  model  to  the  classical  linear  model  is  given  in  Figure  6.  This 
hgure  was  generated  using  A  =  20,  M  =  11  and  N  =  32. 


Figure  6.  Exact  detection  performance  of  the  Classical  Linear  Model  compared 
against  the  Total  Energy  Detector. 

In  situations  where  the  possible  target  location  region  is  quite  large,  we  have 
noticed  that  this  detector  performs  almost  as  well  as  the  GLRT  derived  in  Section 
2.3.  This  is  because  the  total  energy  detector,  unlike  the  GLRT,  does  not  depend 
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on  the  target  location  estimate.  When  the  possible  target  location  region  is  large, 
using  a  similar  explanation  given  in  [11],  it  can  be  shown  that  the  Pfa  increases 
with  the  number  of  “bins”  searched  for  target  location.  This  translates  as  a  dete¬ 
rioration  of  the  performance  of  the  GLRT.  While  the  total  energy  detector  does 
not  have  this  degradation,  it  cannot  be  used  for  target  localization  like  the  GLRT. 
Figure  7  shows  the  comparison  of  the  performance  of  the  GLRT  against  the  total 
energy  detector  where  we  restricted  the  grid  search  for  the  maximum  eigenvalue 
to  only  a  region  of  3kmx3km  square  around  the  actual  target  location. 


ROC  curves  for  average  ENR  =  10dB  per  sensor 


Figure  7.  Gomparison  of  GLRT  against  the  Total  Energy  Detector  for  observations 
collected  at  11  sensors  with  unknown  time  delays  and  Doppler  shifts 


2.6  Conclusion 

We  have  modeled  the  problem  of  detection  of  LPI  signals  and  derived  the 
GLRT  detector  for  the  model.  We  have  shown  that  the  GLRT  uses  both  the 
energy  and  the  cross-sensor  correlation  and  thus  outperforms  the  currently  used 
detectors  which  only  use  either  the  energy  or  the  cross  correlation.  While  Ending 
the  GLRT,  the  MLE  of  the  target  location  and  velocity  are  also  simultaneously 
computed.  We  also  derived  simpler  detectors  by  making  various  assumptions  in 
the  main  model.  We  have  shown  that  when  a  target  has  to  be  detected  in  a  large 
region,  the  total  energy  detector  is  only  slightly  poorer  than  the  GLRT.  In  the  next 
paper  we  will  discuss  the  performance  of  the  MLE  of  the  target  location  mentioned 
here.  We  are  also  investigating  the  possibility  of  extending  this  technique  to  the 
detection  of  multiple  targets. 
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3  Localization 


3.1  Introduction 

Passive  localization  has  been  used  for  many  years  and  has  always  been  an 
important  topic  of  research  [17,  18,  19,  20].  Localization  can  be  performed  using 
one  or  more  of  the  emitter  location  dependent  properties  of  the  signal  such  as 
angle  of  arrival,  TDOA,  FDOA  or  the  energy  of  the  received  signal.  Over  the 
years  the  general  approach  to  localization  using  the  TDOAs,  commonly  referred 
to  as  the  TDOA  technique,  has  been  to  first  estimate  the  difference  in  the  times 
of  arrival  of  the  signal  at  a  particular  pair  of  sensors  and  then  use  these  TDOAs 
to  estimate  the  location  of  the  emitter.  Knapp  and  Carter  [21,  22]  proposed  a 
generalized  correlation  method  for  the  estimation  of  the  time  delay  for  stationary 
and  relative  motion  cases.  They  modeled  the  signal  as  a  stationary  Gaussian  ran¬ 
dom  process.  Stein  [9]  on  the  other  hand  modeled  the  signal  as  deterministic  but 
unknown  and  derived  the  MLE  for  the  differential  delay  and  Doppler  for  a  two 
sensor  case.  Under  similar  assumptions  for  the  signal,  Yeredor  and  Angel  [23]  have 
derived  the  CRLB  for  the  TDOAs.  After  the  TDOAs  are  estimated  they  are  used 
to  estimate  the  location  of  the  source  [24,  25,  26,  27].  Quite  often,  due  to  network 
capacity  and  computational  constraints,  not  all  sensor  pair  combinations  are  used. 
Fowler  [28]  addressed  the  problem  of  optimal  selection  of  a  subset  of  the  sensor 
pairs.  Torrieri  [17]  proposed  a  linear  least  squares  estimator  where  the  nonlinear 
relation  between  the  TDOAs  and  the  emitter  location  is  linearized  by  expanding 
it  in  a  Taylor  series  about  a  reference  point  and  retaining  the  hrst  two  terms.  This 
is  an  iterative  method  which  requires  some  kind  of  a  priori  information  in  order 
to  obtain  an  initial  guess.  Alternatively,  Chan  and  Ho  [18]  use  an  intermediate 
variable,  which  is  a  function  of  the  emitter  location,  in  order  to  linearize  the  non¬ 
linear  equations.  They  use  a  two-step  weighted  least  squares  (WLS)  algorithm. 
Additionally,  when  the  signal  waveform  is  known,  localization  may  be  performed 
from  the  times  of  arrival  (TOAs)  instead  of  the  TDOAs  [29,  30].  In  such  TOA 
based  techniques  the  unknown  transmission  time  occurs  as  a  nuisance  parameter 
which  will  have  to  be  estimated.  Do  et.  al  [31]  have  shown  that  the  TOA  and 
the  TDOA  measurements  are  transformable  to  each  other  without  a  loss  of  in¬ 
formation  regarding  positioning  and  thus  the  position  estimations  based  on  them 
should  be  theoretically  equivalent.  The  above  techniques  may  be  called  two-step 
techniques  because  the  TOAs/TDOAs  are  first  estimated  at  the  local  sensors  and 
these  TOA/TDOA  estimates  are  used  in  a  second  step  to  compute  the  location  of 
the  emitter. 

Weiss  and  Amar  [32,  33,  34,  35]  have  shown  that  the  two-step  approach  is 
sub-optimal  and  proposed  a  direct  position  determination  (DPD)  approach.  Weiss 
had  derived  the  MLE  for  the  source  location  for  the  case  of  a  stationary  narrow- 
band  radio  frequency  transmitter  using  multiple  stationary  receivers  in  [32].  He 
uses  a  continuous  time  model  and  quickly  considers  the  sampled  version  without 
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discussing  the  effects  of  sampling  on  the  emitter  location  estimate.  He  shows  that 
the  MLE  of  the  emitter  location  is  obtained  by  maximizing  a  quadratic  form  of 
the  signal  samples  whose  coefficients  are  functions  of  the  emitter  location.  He 
considers  the  two  cases  of  signal  known  and  signal  unknown  but  leaves  out  a  more 
important  case  -  signal  known  but  transmission  time  unknown  which  is  most  likely 
to  occur  in  real-world  situations.  He  does  not  discuss  the  Cramer-Rao  lower  bound 
(CRLB)  for  this  problem.  There  is  an  inherent  ambiguity  in  the  commonly  used 
model  and  Weiss  uses  a  constraint  on  the  signal  samples  to  resolve  the  ambiguity. 
No  discussion  is  provided  on  the  generality  of  the  constraint  as  to  why  it  is  an 
appropriate  constraint,  how  it  resolves  the  ambiguity  and  whether  it  reduces  the 
performance.  In  [33]  Amar  and  Weiss  extend  the  approach  to  a  multiple  emitters 
case.  In  [34]  they  address  the  problem  of  localization  using  only  the  Doppler  fre¬ 
quency  shifts  and  in  [35]  they  consider  the  case  of  a  single  stationary  emitter  and 
moving  receivers.  In  all  these  cases  the  results  are  similar,  i.e.,  the  MLE  for  the 
emitter  location  is  obtained  by  maximizing  a  quadratic  form  of  the  signal  whose 
coefficients  are  functions  of  the  emitter  location.  The  derivation  of  CRLB  is  at¬ 
tempted  in  the  later  papers  but  is  not  sufficiently  simplihed.  The  effect  of  sampling 
the  signal  is  not  discussed  in  any  of  the  papers.  Similar  constraints  are  used  to 
resolve  the  ambiguity  in  the  following  papers  but  no  discussion  is  provided  on  the 
effects  of  the  constraint. 

In  this  paper  we  consider  the  case  of  a  single  stationary  emitter  and  a  net¬ 
work  of  stationary  receive  sensors.  We  address  many  of  the  short-comings  of 
[32,  33,  34,  35].  We  use  a  continuous  time  model  and  provide  a  straightforward 
derivation  for  the  MLE  of  the  emitter  location  for  the  two  cases  of  signal  wave¬ 
form  known  with  unknown  transmission  time  and  signal  waveform  unknown  with 
unknown  transmission  time.  Our  model  is  valid  for  either  narrow-band  or  broad¬ 
band  signals,  lowpass  or  highpass  signals.  We  discuss  the  effect  of  sampling  the 
signal  on  the  emitter  location  estimate.  Using  simulations,  we  compare  the  MLE 
against  a  conventional  TDOA  technique.  We  show  that  the  variance  of  the  MLE 
is  two  to  three  orders  of  magnitude  lower  than  the  conventional  TDOA  technique. 

A  more  difficult  problem  is  deriving  the  CRLB.  If  the  signal  waveform  is 
assumed  unknown  along  with  the  time  of  arrival  (TOA)  and  the  attenuation  fac¬ 
tor,  then  the  commonly  used  model  has  an  ambiguity.  This  ambiguity  comes  to 
light  when  deriving  the  CRLB.  Because  all  the  unknowns  in  the  model  cannot 
be  uniquely  resolved,  the  Fisher  information  matrix  (FIM)  becomes  singular.  We 
address  this  ambiguity  in  detail  and  derive  the  necessary  steps  to  remove  it.  Then 
we  derive  the  non-singular  FIM.  The  inverse  of  the  FIM  is  the  CRLB.  CRLB  gives 
the  theoretical  lower  bound  on  the  variance  of  any  unbiased  estimator. 

An  important  application  of  the  CRLB  is  in  deriving  an  optimal  sensor  con- 
hguration.  The  performance  of  a  location  estimator  depends  on  the  placement  of 
sensors.  A  particular  conhguration  of  the  sensors  is  called  optimal  if  it  optimizes  a 
norm  of  the  FIM.  A  quite  common  result  [36,  37]  is  to  place  the  sensors  around  the 
emitter  in  an  equi-angular  conhguration.  But  when  the  sensors  are  geographically 
constrained  the  problem  becomes  much  more  difficult.  We  introduce  this  problem 
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in  Section  4.1  and  provide  optimal  sensor  configurations  for  the  three  and  four 
sensor  cases. 

In  Section  6.2  we  provide  a  detailed  description  of  the  problem.  In  Section 
3.3  we  give  the  CRLBs  and  the  MLEs.  Here  we  analyze  the  case  of  signal  wave¬ 
form  unknown  and  the  special  case  of  signal  waveform  known,  both  cases  with 
an  unknown  transmission  time.  In  Section  6.4  we  use  Monte  Carlo  simulations 
to  compare  the  performance  of  the  MLE  against  the  conventionally  used  TDOA 
technique.  We  show  that  at  higher  signal-to-noise  ratios  (SNRs)  the  variance  of  the 
MLE  approaches  the  CRLB.  In  section  4.1  we  introduce  the  problem  of  optimal 
sensor  configurations  and  give  some  results.  Conclusions  are  provided  in  Section 
3.5.  Most  of  the  mathematics  is  provided  in  the  appendices.  In  Appendix  B.l  we 
derive  a  compact  expression  for  the  FIM.  Appendix  B.2  has  the  derivation  of  the 
MLE.  Appendix  B.3  presents  the  properties  of  a  matrix  we  use  in  the  model.  In 
Appendix  B.4  we  discuss  the  transformation  of  the  parameters  and  the  constraints 
used  in  order  to  remove  the  ambiguity  in  the  model.  In  Appendix  C  we  derive 
the  unconstrained  optimal  sensor  configuration  and  the  constrained  optimal  sensor 
configurations  for  the  three  and  four  sensor  cases. 


3.2  Methods,  Assumptions,  and  Procedures 

Suppose  that  a  stationary  emitter  is  located  at  an  unknown  location 
and  a  network  of  M  sensors  are  located  at  known  locations  {xi,  i/i),  i  =  0, 1,  •  •  •  ,  M— 
1  as  shown  in  Figure  27.  For  simplicity  we  are  assuming  a  two  dimensional  case. 
Extension  to  the  three  dimensional  case  is  straightforward.  The  sensors  are  all 
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Figure  8.  Physical  placement  of  the  sensors  (for  M=4)  and  the  Emitter  position 
used  for  simulation. 

synchronized  in  time  and  each  of  the  sensors  intercepts  the  signal  within  the  time 
interval  (0,T).  The  emitter  transmits  an  unknown  signal  s{t)  for  an  unknown 
duration  Tg  <  T  starting  at  an  unknown  time  <  T  .  We  shall  assume  that 
the  transmitted  signal  s{t)  is  real.  It  can  be  narrowband  or  wideband,  lowpass 
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or  bandpass.  After  interception,  the  signal  received  at  sensor  i  in  the  presence  of 
noise  can  be  written  as 


Tiit)  =  Ais{t  —  Ti)  +  Wiit),  0<t<T,  i  =  0, —  1  (17) 

where  Wi{t)  is  a  zero  mean  wide  sense  stationary  additive  white  Gaussian  ran¬ 
dom  process  with  spectral  density  Aj’s  are  the  unknown  attenuations  due  to 
propagation  loss,  assumed  real,  and  the  r^s  are  the  unknown  TOAs  given  by 


A  = 


-  XiY  +  -  ViY 


+  to)  *  —  0)  1)  ■  ■  ■  I  ^  1 


(18) 


where  c  is  the  propagation  speed  of  the  signal.  We  assume  that  the  noise  at  a 

sensor  is  independent  of  the  noise  at  any  other  sensor,  i.e.,  Wi(t)  and  Wj(t)  are 

independent  for  i  j  and  that  the  noise  spectral  density  at  all  the  sensors  is  equal 

to  .  If  the  noise  does  not  satisfy  these  conditions  then  the  problem  becomes 

more  complex.  For  example  if  the  noise  spectral  density  is  different  at  each  sensor 

but  known,  then  the  noise  term  does  not  factor  out  as  in  equation  (22)  but  instead 

exists  in  each  term.  A  more  difficult  problem  is  when  the  noise  spectral  density  is 

different  at  each  sensor  and  unknown,  in  which  case,  the  noise  spectral  densities 

at  each  of  the  sensors  need  to  be  estimated  as  well.  To  keep  the  derivations  simple 

we  assumed  the  above  conditions  for  the  noise.  Notice  that  here  we  do  not  assume 

as  in  [9],  that  r*  <<  T.  Instead  we  just  assume  that  max(ri  —  Tj)  <  (T  —  Tg). 

id 

That  is,  we  are  only  assuming  that  the  observation  interval  is  large  enough  so 
that,  within  the  observation  interval,  the  signal  reaches  both  the  nearest  and  the 
farthest  sensors  from  the  emitter.  Based  on  the  sensor  geometry  it  is  possible  to 
hnd  a  sufficient  condition  on  the  length  of  the  observation  interval.  If  dmax  is  the 
distance  between  the  farthest  pair  of  sensors,  then  the  observation  interval  must 
be  greater  than 

Sampling  the  signal  in  time  has  a  quantization  like  effect  on  the  estimate  of  the 
emitter  location.  This  is  because  if  the  signal  is  sampled,  then  the  TOA  estimates 
are  integer  multiples  of  the  sampling  interval  and  hence  quantized.  For  example 
if  the  signal  is  sampled  at  a  frequency  of  Fg  samples/sec,  then  the  estimate  of 
Ti  is  quantized  with  a  maximum  quantization  error  of  This  can  introduce  a 

maximum  quantization  error  of  ^  in  the  range  estimate.  Therefore,  it  is  possible 
that  the  signal  may  have  to  be  sampled  at  a  rate  much  higher  than  the  Nyquist 
rate  in  order  to  achieve  a  desired  precision  in  the  location  estimate.  We  have  used 
a  continuous  time  model  to  avoid  this  problem  in  our  analysis  and  to  allow  future 
studies  of  errors  due  to  time  synchronization  effects. 

Figure  9  shows  the  signals  received  at  the  four  sensors  shown  in  Figure  27  from 
an  emitter  located  at  (130,  75)km  transmitting  a  Gaussian  chirp.  The  propagation 
loss  was  modeled  as  a  ^  attenuation  in  the  amplitude  of  the  received  signal,  where 
R  is  the  range.  So,  the  farthest  sensor  has  the  largest  TOA  and  smallest  amplitude. 
We  are  assuming  that  the  signal  lies  inside  the  observation  interval  (0,T).  So,  we 
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Signals  received  at  the  sensors  (  no  noise) 


Figure  9.  Signals  received  at  the  four  sensors  when  a  Gaussian  chirp  is  transmitted 
by  the  emitter  located  at  (130,75)km. 


can  assume  that  the  unknown  signal  is  periodic  with  period  T  and  write  it  in  terms 
of  its  Fonrier  series  as 


s{t) 


Qq 

71 


oo 

COS  2'n'nFot  +  bn  sin  27rnFot) 

n=l 


(19) 


where  Fq  =  ^  and  the  Fourier  coefficients  are  given  by 

do  =  X  fo  ~  T  fo  cos27rnFot  dt,  bn  =  Jq  s{t)  sin27mFot  dt 

(20) 

We  are  using  ^  for  the  d.c  component  instead  of  the  standard  Oq  becanse  it 
simplihes  certain  terms  in  the  derivation  of  the  CRLB.  For  a  band-limited  signal 
only  a  hnite  nnmber  of  the  Fourier  coefficients  are  non-zero.  If  the  signal  is  a 
lowpass  signal,  there  exists  an  integer  N  such  that  the  Fourier  coefficients  are  all 
zero  for  n  >  N  and  if  the  signal  is  a  bandpass  signal,  then  there  exist  integers  Ni 
and  iVg,  Ni  <  N2,  snch  that  the  Fonrier  coefficients  are  zero  for  n  <  Ni  and  for 
n  >  N2.  So  we  can  approximate  the  lowpass  signal  s(t)  as  (  for  a  bandpass  signal 
the  snmmation  is  from  Ni  to  iVg  ) 


s{t) 


Qq 

71 


N-l 

(a„  cos  27rnFo^  +  bn  sin  2'KnFQt) 

n=l 
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This  is  an  important  step  as  it  allows  us  to  model  any  unknown  signal  and  reduce 
it  to  a  parameter  estimation  problem.  Now,  if  we  let 

<p  =  [ao  ai  ■  ■  ■  aisf-i  &2  ■  ■  ■  be  the  2iV  —  1  x  1  vector  of  Fourier  coefficients 

and 


h{t) 


■  1 
.71 


cos27rFot  ■■■  cos27r(iV  —  l)Fot 


sin  27rFf)t 


T 


sin  27r(iV  —  l)Fot 


then  we  have  s{t)  =  h^{t)cl).  This  reduces  the  uncountable  unknown  parameter 
set  {s{t)  :  t  G  (0,T)}  to  a  hnite  countable  number  of  unknown  parameters  0. 
Therefore,  we  can  rewrite  the  model  in  equation  (45)  as 


ri{t)  =  Aih^  (t  —  rj)0  +  Wi{t),  0  <t  <T,  i  =  0,l,---,M  —  1  (21) 


Let  T  =  [tq  Ti  •  •  •  and  A  =  [Aq  Ai  ■  ■  ■  Let  0  =  [r'^  4^'^]'^  be 

the  (2M  +  2N  —  1)  x  1  vector  of  unknown  parameters.  If  we  let  rj  =  [x.^  y^,  to]'^ 
and  CK  =  [t]'^  A^  0^]^  then,  using  equation  (46),  we  can  write  the  TOA  vector 
as  a  function  of  ry  as  r  =  g{r]).  So,  the  problem  can  be  stated  as,  given  the 
observations  ri{t),  i  =  -  ,M  —  1  estimate  the  vector  ry.  We  are  only 

interested  in  the  parameters  {Xj,,y.j,)  and  the  rest  of  the  unknown  parameters  are 
nuisance  parameters. 


3.3  CRLB  and  MLE  of  the  Emitter  Location 

Here  we  will  derive  the  CRLB  and  the  MLE  of  the  emitter  location  for  the 
two  cases  of  signal  waveform  unknown  with  unknown  transmission  time  and  signal 
waveform  known  with  unknown  transmission  time. 


3.3.1  Signal  unknown  with  unknown  transmission  time 

For  the  continuous  time  model  in  equation  (45)  the  log-likelihood  function 
[38]  for  sensor  i  is  given  by 


/  =  -—  /  {ri{t)  -  Ais{t  -  Ti)f  dt  (22) 

^^0  Jo 

Since  the  noise  at  different  sensors  is  independent,  and  using  equation  (21)  we  can 
write  the  joint  log-likelihood  function  as 

l{9)  =  /  'Y]  {ri{t)  -  Aih^{t  -  Ti)cl)f  dt  (23) 

^oJo 


The  (2M-I-2A  — 1  x2M-|-2A  — 1)  FIM  [39]  for  this  model  is  given  by  (see  Appendix 
B.1.1) 


(T/2) 

(JVo/2) 


(27rFo)2<^'^LL^<#)(diag(A))2  (27rFo)(<#>^L,^)(diag(A))  (27rFo) (A  ©  A) ,^'^L 

(27rFo)(<#>^L^<^)(diag(A))  {cfiX  A<#)^ 

(27rFo)L'^<#>(A  0  A)^  0A^  (A^A)Im 


(24) 
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where  ©  represents  the  element  by  element  product  (Hadamard  product),  diag(A) 
is  an  M  X  M  diagonal  matrix  with  iih.  diagonal  element  as  Aj,  is  the  M  x  M 
identity  matrix  and  the  {2N  —  1)  x  {2N  —  1)  matrix  L  is  given  by 

[  diag(l,2,...  ,A-1)  J 

-  “[  ^(Af-1,1)  diag(l,2,  •••  ,N  —  1)  ] 

This  FIM  must  be  inverted  in  order  to  hnd  the  CRLB  for  the  unknown  parameter 
vector  0.  By  the  form  of  the  matrix  in  equation  (24)  it  is  easily  shown  (  see 
Appendix  B.4  )  that  the  matrix  is  singular  with  rank  equal  to  two  less  than  full 
rank.  Weiss  [32]  uses  an  ad  hoc  method  to  overcome  this.  We,  however  use 
the  exact  transformation  of  the  parameters  [40]  that  is  required  to  eliminate  the 
singularity  of  the  information  matrix.  The  singularity  arises  because  there  is  an 
ambiguity  in  the  model.  It  is  not  possible  to  uniquely  determine  all  the  unknown 
parameters  in  the  model  in  equation  (45).  This  is  because  of  the  relationship 
between  the  TOA,  attenuation  factor  and  the  signal  waveform.  Suppose  that  in 
equation  (45),  A*  and  s{t)  are  the  true  values  of  the  gain  and  the  signal  waveform 
that  generate  ri{t).  Then  the  pair  of  values  {Ai/'y,  'ys{t))  for  any  non-zero  constant 
7  also  generate  the  same  rj(t).  So,  from  the  observation  rj(t)  it  is  impossible  to 
determine  the  true  values  of  Aj  and  s{t).  A  similar  relationship  exists  between  the 
TOA  and  the  signal  waveform.  Suppose  that  in  equation  (45),  fj  and  s(t)  are  the 
true  values  of  the  TOA  and  the  signal  waveform  that  generate  rj(t).  Then  the  pair 
of  values  ((r^  —  7),  s(t  —  7))  for  any  constant  7  also  generate  the  same  rj(t).  This 
is  more  clearly  demonstrated  in  Figure  10.  Notice  that  given  the  received  signal 
it  is  not  possible  to  determine  whether  the  source  signal  is  Si{t)  with  TOA  ti  or 
S2{t)  with  TOA  T2.  This  causes  the  information  matrix  to  be  at  least  rank  two 


Transmitted  Signal  Received  Signai 


Figure  10.  Ambiguity  when  signal  and  TOA  are  both  unknown. 

dehcient,  as  shown  in  Appendix  B.4.  The  over  parameterization  can  be  resolved 
by  applying  an  appropriate  transformation  that  satishes  certain  constraints  [40]. 
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As  shown  in  Appendix  B.4,  the  appropriate  transformation  for  this  model  is 


T 


/ 


A' 


0' 


=  [  (b  -  To)  {t2  -  To)  ■■■  {tm-1  - 

=  (l/^o)[^l  ■■■ 


1 

0(1,2V-2) 

0(27V-2,1) 

f-V-l 

lAf-l 

Iv-1 

— Itv-i 

diag(ft(-To))0 


(25) 


and  can  he  shown  to  be  the  least  restrictive  constraint  for  identifiahility.  Here 
we  are  using  the  (*)'  notation  to  represent  the  new  parameters  resulting  from 
the  transformation.  Notice  that  the  transformed  parameter  vectors  r'  and  A' 
are  TDOA  and  relative  gain  factor  with  respect  to  sensor  0  and  are  each  reduced 
by  one  parameter  from  r  and  A  respectively  while  the  (p'  is  simply  the  Fourier 
coefficients  of  roif).  Using  these  transformed  parameters  the  model  in  equation 
(21)  can  be  rewritten  as 


roit)  =  {t)4>'  +  Wo{t)  0  <t  <T 

ri(t)  =  A'Ji^ (t  —  t[)4>'  +  Wi{t)^  0  <  t  <  T,  i  =  1,  2,  •  •  •  ,  M  —  1 


where  A'  =  ^  and  t[  =  r— tq.  So,  the  effect  of  the  transformation  is  that  the  signal 
at  sensor  0  is  made  the  reference  signal  and  the  signals  at  all  the  other  sensors  are 
modeled  relative  to  this  reference  signal.  Although  (26)  seems  intuitively  obvious, 
by  arriving  at  it  from  21)  using  the  transformation  in  (25),  we  have  mathematically 
verihed  that  (26)  is  indeed  the  correct  model  to  use  for  the  problem  of  localization 
under  the  unknown  signal  case.  Weiss  [32,  33,  34,  35]  uses  (26)  directly  without 
this  rigorous  argument.  This  is  a  subtle  but  important  result  which  is  overlooked 
by  Weiss.  Now,  let  O'  =  [t''^  A'^  the  (2M  +  2N  —  3)  x  1  vector  of 

the  unknown  transformed  parameters,  r]'  =  [x^,y^Y  ■ 

Notice  that  the  TDOA  vector  r'  is  a  function  of  only  {x^,yY,  he.  r'  = 

The  unknown  transmission  time  to  does  not  appear  and  thus  is  not  a  nuisance 
parameter.  The  problem  is  now,  given  the  observations  rpt),  i  =  0,1,  ■  ■  ■  ,  M  — 
1  estimate  the  vector  r}'.  The  log-likelihood  function  for  this  model  with  the 
transformed  parameters  is  given  by 


*(®')  =  “A  F  ('■»(*)  -  Io  E"i  '  ('■((«)  -  -  T04>'f  dt 

(27) 

As  shown  in  the  Appendix  B.1.1,  the  FIM  for  this  transformed  parameter  vector 

is 

XY  =  HXJh^  (28) 

where  H  is  given  in  (B.29),  or  equivalently. 


T/2 

^  (JVo/2) 


(2tvFo)^cI>'  ^LL7’,),'(diag(A'))^  (27rFo)(<#>'  ■^L<#.')(diag(A')) 

(27rFo)(0'  '^L7’<^')(diag(A'))  (</>'  ^<)>')I(M-1) 

(27rFo)L7’,),'(A' 0  A')^  <#>'A' 7’ 


(2,rFo)(A' 0  A')<#>' '^L  ' 

A'<#>'  ^ 


(29) 
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The  FIM  for  the  corresponding  vector  ol'  is  given  by  [39] 


X„/  = 


do'  V  ^  f  d9'  \ 

dcx'T)  '^'\dcx'T) 


/  de' 


(HijHy-' 


dG' 

dct''^ 


(30) 


(do'  \ 

— — yj;  j  is  given  in  equation  (B.18).  In  Appendix  B.2.1 

we  show  that,  the  MLE  for  the  emitter  location  is  obtained  by  maximizing  over 
the  maximum  eigenvalue  of  the  M  x  M  cross-correlation  matrix  B'  = 

Y-y'T  =  E£o'  yxj  where  Y'  =  [  /„  y)  ...  yV_,  ]  with 

y'o  ^  Mi)h(t)  dt  and  y)  = rj(t)h(i  -  t')  dt,  i  =  1, 2, . . .  M  -  1.  (31) 


That  is, 
or  equivalently. 


f)'  =  argmaxAmax(B') 

rji 


(32) 


(dr,  yr)  ^  arg  max  An,aa(B')  (33) 

(xT,yT) 

where  Amax  represents  the  maximum  eigenvalue.  The  matrix  B'  is  real  symmetric 
and  positive  dehnite  and  so  the  maximum  eigenvalue  is  real  and  positive. 


3.3.2  Signal  known  with  nnknown  transmission  time 

Quite  often  in  practical  situations  it  is  possible  that  the  signal  waveform  is 
known  but  the  exact  transmission  time  to  is  unknown.  In  this  case  the  number  of 
unknowns  is  reduced  to  2M.  Let  ^  =  [r^  be  the  2M  x  1  unknown  parameter 
vector.  Similar  to  (23),  the  log-likelihood  function  is  given  by 

^  M-l 

/(C)  =  “  Aih^{t  -  dt  (34) 

/Vo  Jo 


where  cf)  is  known.  The  FIM  for  this  model  is  given  by  (see  Appendix  B.1.2) 


.  (r/2) 


(2irF(i)^</)’’LL’’(/)(diag(A))^  (2irFo)(i))’’L</))(diag(A)) 
(2irF„)(d>’’L’'d>)(diag(A))  (0^<A)Im 


(35) 


This  matrix  is  not  singular  because,  for  this  case,  the  unknown  parameters  in 
the  model  can  be  uniquely  determined.  Therefore,  there  is  no  need  to  transform 
the  parameters  as  in  the  case  of  the  unknown  signal.  Although,  the  unknown 
transmission  time  to  is  still  retained  here  as  the  nuisance  parameter.  For  this 
model  it  is  shown  in  Appendix  B.2.2  that  the  MLE  for  emitter  location  and  the 
unknown  transmission  time  is  given  by 

^  =  argmax  (/)^B(/)  (36) 

G 
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where  B  =  YY^  and  Y  =  [yg  ■  ■  ■  y^-i]  with  y^  =  ri{t)h{t  -  Ti)  dt,  i  = 
0, 1,  •  •  •  M  —  1.  B  is  a  function  of  {x^,  to).  Using  the  fact  that  {t  —  Ti)(t)  = 
s{t  —  Ti)  and  T]  =  [xj,  y^,  we  can  rewrite  equation  (36)  as 

M-l  /  rT  \  2 

{x^,y^,io)  =  argmax  ^  (  /  ri{t)s{t  -  n)  dt  j  (37) 

Equation  (37)  is  simply  the  correlation  values  between  the  known  signal  waveform 
and  the  observed  signal  at  each  of  the  sensors,  summed  over  all  the  sensors.  The 
emitter  location  that  yields  the  values  of  the  TOAs  that  maximize  the  expression 
in  equation  (37)  is  the  MLE  of  the  emitter  location. 


3.4  Results  and  Discussion 

In  order  to  evaluate  the  performance  of  the  MLE  we  have  run  some  simulations 
and  compared  the  performance  against  the  CRLB  and  against  a  typically  used 
TDOA  approach.  The  TDOA  approach  was  implemented  as  a  two-step  algorithm 
where,  in  the  first  step,  the  TDOA  AA  estimates  were  obtained  by  cross-correlating 
the  signal  at  each  sensor  with  the  signal  at  sensor  0.  Then  a  1km  x  1km  region 
around  the  true  emitter  location  was  split  into  100  x  100  grid  points  and  for  each 
emitter  location  on  the  grid  point  the  TDOAs  were  computed  using  the  formula 

_  V  +  {Vt  -  Vif  V  +  {Vt  -  VoY  i  =  l2---  M-l 

*  C  C  1  1  1  ■: 

Next  the  least  squared  error  (LSE)  between  the  estimated  TDOAs  and  the  com¬ 
puted  TDOAs  was  calculated  as 

M-l 

LSE=  5^(AA- Ar,)2 

i=l 


This  LSE  is  a  function  of  the  emitter  location  {x^,y^).  The  emitter  location 
that  minimized  the  LSE  is  the  estimate  of  the  emitter  location.  Since,  we  know 
that  the  LSE  is  a  2-dimensional  parabolic  function  of  the  emitter  location,  we 
improved  the  accuracy  by  htting  a  parabola  through  the  10000  points  at  which 
the  LSE  was  computed.  Then  by  using  the  analytical  formula  for  the  minimum 
location  of  a  2-dimensional  parabola,  we  computed  the  minimum.  To  measure  the 
levels  of  the  zero  mean  additive  white  Gaussian  noise,  we  used  a  metric  called 
the  average  SNR  (ASNR).  The  ASNR  is  the  ratio  of  the  average  signal  power 
to  the  noise  power  at  each  sensor  averaged  over  all  the  sensors,  i.e.,  if  Vsi  = 
|AjpL  /o  k(^)P  dt  is  the  average  power  of  the  signal  at  the  fth  sensor  and  ^  is 
the  noise  spectral  density  at  the  fth  sensor,  then  the  SNR  averaged  over  M  sensors 

is  given  by  10 log  {No/2)If  /2) )  11  (a)  shows  a  realization 

of  the  LSE  function.  For  the  simulation  we  have  used  4  sensors  placed  at  the 
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y(km)  74.5  129.5  x  (km) 
x10“^ 


y(km)  74.5  129.5  x  (km) 


Figure  11.  (a)  A  realization  of  the  LSE  function  at  an  ASNR  =  —10  dB.  (b)  A 
realization  of  the  Likelihood  function  at  an  ASNR  =  —10  dB. 


coordinates  shown  in  Figure  27  and  the  emitter  was  placed  at  the  coordinates 
(130,  75)  km,  also  as  shown  in  Figure  27.  The  integrations  in  (20)  and  (31)  are 
approximated  using  summations  with  St  =  0.33  ns  which  means  the  sampling 
frequency  is  Fg  =  300  MHz.  The  reason  for  choosing  such  high  sampling  frequency 
is  that  at  this  frequency,  the  position  quantization  error  due  to  sampling  is  in  the 
order  of  {c/Fg)  =  10“^  km.  A  Gaussian  chirp  dehned  by 


s(t)  =  exp 


sin(27rmt^) 


was  used  as  the  unknown  transmitted  signal  waveform.  Figure  12  shows  the  trans¬ 
mitted  signal  waveform.  Notice  that  the  signal  is  assumed  to  be  approximately 
zero  for  t  <  0  and  for  t  >  Tg.  We  set  Tg  =  5/is  and  =  200,  OOOtt.  The  obser¬ 
vation  interval  at  each  of  the  sensors  was  taken  to  be  T  =  0.2  ms.  The  unknown 
transmission  time  of  the  signal  was  set  to  to  =  0.07  ms.  With  this  conhguration  the 
maximum  TDOA  is  0.0988  ms.  The  frequency  spectrum  of  the  Gaussian  window 
is  given  by  |5'(F)|  =  (\/^/o'^)exp  (— 27r^F^/(T^)  [39],  and  the  bandwidth  of  the 
chirp  is  BW  =  1.5  MHz.  The  rate  of  change  of  frequency  for  the  linear  chirp  was 
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Transmitted  Signal 


Figure  12.  Transmitted  signal  waveform. 


chosen  to  be  m 


BW 

T 

s 


3  X  10®.  A  plot  of  the  Fourier  coefficients  of  the  signal  is 


shown  in  Figure  13.  A  value  of  A^  =  2  x  BW  xT  =  600  was  used  to  have  a  total  of 
2N  —  1  =  1199  unknown  Fourier  coefficients.  Notice  that  the  Fourier  coefficients 
are  almost  zero  for  n  >  600.  We  set  the  ASNR  at  —20  dB  and  ran  a  total  of  300 


Fourier  Coefficients 


Figure  13.  Fourier  coefficients  plot. 

Monte  Carlo  simulations  to  generate  the  scatter  plot  and  the  corresponding  95% 
error  ellipse  which  are  shown  in  Figure  14.  This  is  also  called  the  95%  conhdence 
ellipse.  That  is,  if  this  estimator  is  used  a  large  number  of  times  for  localization, 
then  around  95%  of  those  times  the  true  location  of  the  emitter  will  lie  within  this 
ellipse.  To  compute  the  MLE  the  maximization  was  performed  using  a  grid  search. 
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Figure  14.  Scatter  plot  and  the  corresponding  95%  error  ellipse  of  the  MLE  for  an 
ASNR  =  -20  dB. 

We  used  a  grid  of  size  1  km  x  1  km  with  the  grid  points  0.01  km  apart  to  have 
a  total  of  100  x  100  =  10000  points.  The  grid  is  shown  with  a  dotted  line  in  the 
hgure.  At  —20  dB,  the  variances  of  the  MLEs  of  were  (0.0021,  0.0006)  km^ 

and  the  respective  CRLBs  were  (0.0012,  0.0003)  km^.  The  hgure  does  not  show  300 
points  because  some  points  lie  on  top  of  the  others  due  to  position  quantization 
induced  by  the  hnite  number  of  grid  points  in  the  grid  search  for  maximization. 
Due  to  the  complex  nature  of  the  likelihood  function  (see  Eigure  11  (b)),  it  is  not 
possible  to  use  any  curve  htting  techniques  to  reduce  this  quantization  effect  as  in 
the  case  of  the  conventional  TDOA  approach. 

Figure  15  shows  the  comparison  of  the  variances  of  the  MLE  and  the  typical 
TDOA  approach  against  the  CRLB  for  different  SNR  values.  Notice  that  for  the 
ASNR  values  below  —30  dB,  the  variance  of  the  MLE  remains  hat.  This  is  because 
of  the  restriction  imposed  by  the  hnite  grid  size.  As  the  ASNR  increases  above 
—30  dB,  the  variance  of  the  MLE  reduces  rapidly  to  approach  the  CRLB  at  around 
—  10  dB.  Due  to  the  nature  of  the  conventional  TDOA  approach  it  breaks  down 
for  the  ASNR  values  below  —17  dB.  So  this  hgure  has  the  variances  of  the  TDOA 
approach  only  for  the  average  SNR  values  above  —17  dB.  On  the  other  hand,  for 
this  particnlar  setnp,  the  resnlts  for  the  MLE  are  reliable  for  the  ASNR  valnes  as 
low  as  —30  dB.  It  is  qnite  obvious  from  this  hgure  that  performance  of  the  MLE  is 
very  mnch  better  than  a  typical  TDOA  approach.  In  this  case,  the  MLE  performs 
as  good  as  a  typical  TDOA  approach  for  an  ASNR  value  of  about  10  dB  less  than 
that  for  the  TDOA  approach.  Also  notice  that  at  around  —10  dB,  the  variance  of 
the  MLE  is  almost  two  orders  of  magnitnde  less  than  that  of  the  TDOA  approach. 
We  have  noticed  that  for  certain  sensor-emitter  conhgnrations,  particularly  when 
the  emitter  was  very  close  to  the  sensors,  the  variance  of  the  MLE  is  np  to  three 
orders  of  magnitnde  less  than  that  of  the  TDOA  approach.  Therefore,  nnder  LPI 
scenarios  where  a  conventional  TDOA  techniqne  cannot  be  reliably  used,  the  MLE 
can  be  used. 
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CRLB  vs  MLE  and  TDOA  variance  for  different  SNR 


-35  -30  -25  -20  -15  -10 

SNR  (dB) 

Figure  15.  Comparison  of  variances  for  the  emitter  location  estimate  using  MLE 
and  a  typical  TDOA  approach  against  the  CRLB  for  different  average  SNR  values. 

3.5  Conclusions 

We  have  derived  a  direct  positioning  estimator  for  an  emitter  location.  This 
is  the  maximum  likelihood  estimator.  We  have  shown  that  for  an  unknown  signal 
case,  the  model  that  is  conventionally  used  has  an  inherent  ambiguity  and  so  all  the 
unknown  parameters  cannot  be  uniquely  determined.  We  derived  an  appropriate 
transformation  of  the  parameters  and  re-parameterized  the  model  to  remove  the 
ambiguity.  We  have  shown  that  for  the  special  case  of  a  known  signal  with  unknown 
transmission  time,  there  is  no  ambiguity  in  the  model.  We  derived  the  MLE  and 
the  FIM  for  the  model.  The  performance  of  the  MLE  was  compared  against  a 
typical  two-step  TDOA  based  localizer  and  against  the  CRLB.  The  performance 
of  the  MLE  is  signihcantly  better  than  a  typical  two-step  TDOA  based  localizer. 
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4  Optimal  Sensor  Configuration 


4.1  Introduction 

Apart  from  helping  evaluate  the  performance  of  the  MLE,  the  CRLB  has  other 
important  applications.  One  such  application  is  to  determine  the  optimal  sensor 
configurations.  From  equations  (30)  and  (B.18)  we  see  that  the  FIM  of  the  emitter 

dr' 


location  vector  [x^  depends  on  the  Jacobian 


dr)‘ 


J  T 


Now,  this  Jacobian 


is  dependent  on  the  sensor  emitter  geometry.  Figure  16  shows  different  geometries 
and  the  corresponding  CRLBs  for  an  ASNR  of  —20  dB.  It  can  be  noticed  that  for 
configurations  that  surround  the  emitter,  the  CRLB  is  smaller.  The  problem  of 


Geometry  1  Geometry  2 


Geometry  3  Geometry  4 


Figure  16.  Different  Sensor  Configurations. 


optimizing  the  sensor  configuration  falls  under  a  broader  class  of  problems  called 
“Optimal  Design  of  Experiments” .  Generally,  a  norm  of  the  FIM  is  chosen  as  the 
optimization  criterion  and  depending  on  the  norm,  the  corresponding  configura- 


Table  2.  CRLB  for  different  geometries 


Geometry 

GRLB(x) 

GRLB(y) 

1 

0.0012 

0.0006 

2 

0.0024 

0.0021 

3 

0.0015 

0.0001 

4 

0.0656x10-^ 

0.2852x10-^ 
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tion  is  classified  as  A-optimal,  where  the  trace  of  the  inverse  of  the  information 
matrix  is  minimized,  D-optimal,  where  the  determinant  of  the  information  matrix 
is  maximized,  T-optimal,  where  the  trace  of  the  information  matrix  is  maximized 
etc.  Lee  [36]  has  derived  the  A-optimal  sensor  conhgurations  for  the  conventional 
TDOA  approach.  He  has  shown  that  the  optimal  conhguration  is  to  place  the  sen¬ 
sors  with  an  equi-angular  spacing  around  the  emitter.  Similar  results  were  derived 
by  Bishop  et  al  [37]  for  the  D-optimal  conhguration.  We  can  notice  consistent 
results  in  the  example  in  Figure  16.  When  the  sensors  are  around  the  emitter  as 
in  Case  4,  the  CRLB  is  greatly  reduced.  Placing  the  sensors  around  the  emitter 
may  not  however  be  practically  feasible.  A  practical  problem  of  interest  is  to  hnd 
an  optimal  sensor  conhguration  where  the  sensors  are  restricted  to  a  sector.  The 
A-optimal  conhguration  for  this  problem  is  addressed  by  Lee  [36]. 


4.2  Methods,  Assumptions,  and  Procedures 

Now,  from  (C.l),  we  have  the  FIM  for  the  emitter  location  as 


Xri'  — 


dr' 
dr)'  ^ 


(27rFo)>''^'LL'^>'(diag(A'))" 


dr' 
dr)'  ^ 


If  we  assume  that  the  signal  level  is  the  same  at  all  sensors,  (Aj  =  A)  for  all  i,  we 
have  the  above  FIM  as 


=  (27rFo)>'^LL">'^ 


dr' 
dr)'  ^ 


dr' 

dr)''^ 

dr' 


The  geometry  dependent  term  in  the  above  equation  is  (  ^ '  „  )  ( 

\dr)'  ^  J  \dr)'  ^ 

fore,  an  emitter-sensor  conhguration  that  maximizes  the  determinant  of 


dr' 


dr)- 


,1  T 


dr' 
dr)''^ 


is  called  the  D-optimal  conhguration. 


There- 


We  dehne  the  unconstrained  problem  as  one  where  the  sensors  are  allowed  to 
be  anywhere  in  angle  (  equivalent  to  being  on  a  circle  )  and,  in  the  constrained 
conhguration  problem,  the  sensors  can  only  he  on  an  arc  of  the  circle. 


4.3  Results  and  Discussion 

We  have  derived  the  D-optimal  conhguration  for  three  and  four  sensor  cases 
(see  Appendix  C).  That  is,  we  derived  an  emitter-sensor  conhguration  that  max¬ 
imizes  the  determinant  of  the  FIM  of  the  emitter  location.  In  Appendix  C.l  we 
show  that  for  the  unconstrained  problem,  a  D-optimal  conhguration  is  to  place  the 
sensors  at  equi-angular  distance  around  the  emitter.  This  appears  rather  contra¬ 
dictory  since  it  is  required  to  know  the  location  of  the  emitter  in  order  to  optimally 
conhgure  the  sensors  for  localization.  In  fact,  this  result  has  applications  in  mo¬ 
bile  emitter  localization  where  in  the  recursive  tracking  of  the  emitter  using  a 
sequence  of  TDOA  measurements  is  performed  [41].  This  result  is  consistent  with 
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the  results  of  Lee  and  Bishop.  Our  proof  in  Appendix  C.l  is  much  simpler  than 
previously  published  results.  But  this  conhguration  is  not  the  unique  D-optimal 
conhguration.  As  the  number  of  sensors  increases,  there  are  many  more  conhg- 
urations  which  are  D-optimal  but  are  difficult  to  hud.  For  example,  in  the  case 


TT 


of  six  sensors,  placing  the  sensors  at  an  angular  distance  of  —  radians  from  each 


other  on  the  circle  around  the  emitter,  as  shown  in  Figure  17  (a),  is  one  D-optimal 
conhguration.  Another  D-optimal  conhguration  is  to  place  two  sensors  at  each 
of  the  vertices  of  an  equilateral  triangle  that  is  inscribed  in  the  circle  around  the 
emitter,  as  shown  in  Figure  17  (b).  Here  the  sensors  are  supposed  to  he  on  top  of 
each  other  but  for  visual  clarity,  they  are  slightly  displaced. 


Figure  17.  Two  D-optimal  conhgurations  for  the  case  of  M=6  sensors. 

The  constrained  problem  is  much  more  difficult  to  solve.  The  setup  of  a 
constrained  problem  is  as  shown  in  Figure  18.  The  sensors  can  lie  only  on  the  arc 
with  half  angle  6  which  is  marked  in  bold  in  the  hgure.  Using  the  transformation 


Figure  18.  Constrained  sensor  conhguration  setup.  The  sensors  can  lie  only  on  the 
arc  marked  in  bold. 
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of  variables  given  in  Appendix  C  we  have  the  determinant  as 


where 


'  dr'  Y 

( ) 

=  det(JJ^) 

dr]'  ^  ) 

\  d-q'  ) 

cos  'Ipo 

cos'll  ' 

•••  COS^/>M-l 

sin'^0 

sin  ipi  ' 

•••  sill'd  M-l 

1 

1 

1 

where  tjji  is  the  angle  of  sensor  i  from  the  positive  x-axis  measured  at  the  emitter 
and  —9  <  ipi  <  6  for  alH  =  0,  •  •  ■  ,  M  —  1.  Notice  that  when  M  =  3,  we  have 


J 


COS'^O 
sin  -ipo 
1 


cos'll 

sin'll 

1 


cos '^2 

sin  1^2 
1 


and  det(JJ^)  is  proportional  to  the  square  of  the  area  of  the  triangle  with  vertices 
at  the  sensor  locations.  So,  in  geometrical  terms  the  problem  can  be  viewed  as 
hnding  the  triangle  with  the  maximum  area  that  has  vertices  on  the  constraining 
arc.  In  order  to  increase  the  area,  the  base  and  the  height  of  the  triangle  must  be 
increased.  So,  intuitively  it  can  be  seen  that  the  optimal  conhguration  for  the  three 
sensor  case  is  as  shown  in  Figure  19  (  A  rigorous  proof  is  provided  in  Appendix 
C.2  ). 


Figure  19.  The  D-optimal  conhguration  for  three  sensors.  The  triangle  with  the 
maximum  area  is  also  shown  here. 

This  geometrical  interpretation  can  be  extended  to  any  number  of  sensors 
using  the  Cauchy-Binet  formula.  For  an  arbitrary  M,  we  have 

det(JJ^)  =  det(J5J^) 

where  [M]  is  the  set  {0, 1, ...,  M  —  1},  and  is  the  set  of  3-combinations  of  [M] 
(i.e.,  subsets  of  size  3).  Js  is  the  3x3  matrix  whose  columns  are  the  columns 
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of  J  at  indices  from  S.  Now,  each  det(J5J^)  is  the  square  of  the  area  of  the 
triangle  formed  by  the  sensors  with  indices  in  S.  So,  maximizing  the  determinant 
is  equivalent  to  maximizing  the  sum  of  the  squares  of  the  areas  of  the  triangles 
formed  by  all  combinations  of  sensors  taken  three  at  a  time.  There  is  a  total  of 
such  combinations.  We  hypothesize  that  in  general  a  D-optimal  conhguration 
for  any  number  of  sensors  will  consist  of  distributing  the  sensors  in  appropriate 
proportions  at  the  center  and  at  the  two  end  points  of  the  arc.  There  are  three 
D-optimal  conhgurations  for  the  four  sensor  case,  which  are  given  in  Figure  20 
(  proof  in  Appendix  C.2).  These  derivations  are  very  specihc  to  the  particular 
number  of  sensors  and  are  hard  to  extend  to  higher  number  of  sensors.  In  a  future 
paper  we  will  discuss  this  constrained  optimal  conhguration  in  detail. 


Figure  20.  The  three  D-optimal  conhgurations  for  four  sensor  case. 


4.4  Conclusions 

We  have  introduced  the  problem  of  optimizing  the  sensor-emitter  geometry. 
We  derived  the  D-optimal  conhguration  for  the  unconstrained  geometry  problem. 
For  the  constrained  geometry  problem,  we  have  derived  the  D-optimal  conhgura¬ 
tions  for  the  three  and  four  sensor  cases. 
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5  Knowledge  Aided  Design 

5.1  Introduction 

In  knowledge  aided  design,  the  knowledge  of  the  terrain  is  used  for  localiza¬ 
tion.  Presence  of  objects  in  a  terrain  induces  some  azimuth  modulation  due  to 
obstruction,  reflection,  diffraction  etc  of  the  signal.  Knowledge  of  the  location  of 
these  obstacles  and  the  type  of  modulation  they  induce,  can  aid  in  the  localization 
of  the  target  of  interest.  This  is  particularly  useful  in  localization  in  the  urban 
environments  where  there  is  buildings,  trees  etc  which  act  as  the  obstacles.  In  this 
section  we  will  derive  the  Fisher  information  matrix  and  show  how  the  azimuth 
modulation  induced  by  an  obstacle  can  increase  the  information. 

5.2  Methods,  Assumptions,  and  Procedures 

Suppose  that  we  have  M  sensors  located  at  {xi,yi)  i  =  1,  •  •  •  ,  M.  Let  an 
emitter  located  at  an  unknown  location  transmit  a  known  signal  waveform 

s{t).  Assuming  that  the  sensors  are  in  a  direct  line  of  sight  from  the  emitter  and 
that  there  is  no  multipath,  the  signal  received  at  the  sensors  is  modeled  as 

ri{t)  =  Ais{t  -  Ti)  +  Wi{t)  i  =  (38) 

where  AjS  are  the  unknown  attenuation  factors  and  r^s  are  the  unknown  time- 
delays.  Wi(t)  is  white  Gaussian  noise.  We  want  to  determine  the  amount  of 
information  in  the  azimuth  modulation. 


5.3  Results  and  Discussion 

The  FIM  for  the  model  in  (38)  is  given  by  (see  Appendix  D  ) 


I'iXj, ,  y^, ) 


_p2 


(Wo/2)  c? 


cos^  'ipi  sin  -ipi  cos  'ipi 
sin  ipi  cos  ipi  sin^  'ipi 


+ 


£ 

(No/2)  2^i=l 


(39) 


Here  the  A^s  and  r^s  are  known  functions  of  the  target  location.  For  example, 
assuming  there  are  no  obstacles,  we  have 


Ai{.x^,y^) 


_ G, _ 

{x^  -  XiY  +  {y^  -  yiY 


Ri 


where  Gi  is  a  constant  and 


=  (l/c)\/(a^T  -  +  {Vt  -  ViY  =  (l/c)i?i 
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where  c  is  the  propagation  speed  of  the  signal.  So,  the  FIM  for  this  problem  is 
given  by  (  see  Appendix  D.l) 


£ 


M 

E 


rx2 

COS^  fji 

sin  -ipi  cos  fji 

'R} 

sin  -ipi  cos 

hi  sin^  Ipi 

+ 


cos^  -ipi  sin  'ipi  cos 
sin  cos  sin^  ipi 


(40) 


—  Xt  + 

where  the  hrst  term  is  the  contribution  from  the  time-delays  and  the  second  term 
is  the  contribution  from  azimuth  modulation.  Note  that  the  information  in  the 
azimuth  modulation  is  greater  than  the  information  in  the  time-delays  if 

F^Gl  & 

That  is,  for  short  ranges  as  in  urban  environments,  or  for  signals  with  narrow 
band-width,  the  information  in  the  azimuth  modulation  beeomes  comparable  to  the 
information  in  the  time- delays.  A  possible  range  and  bandwidth  is  given  in  the 
following  example. 


5.3.1  Example  1  -  No  Obstacles 

We  have  4  sensors  and  an  emitter  located  at  (945,  810)m  as  shown  in  Figure 
21.  The  transmitted  signal  is  a  Gaussian  pulse  given  by 


s{t)  =  exp 


The  length  of  the  signal  in  time  is  Tg  =  50/is  and  its  root-mean-square  (RMS) 
bandwidth  is  =  111.07  kHz.  This  is  also  shown  in  Figure  21.  So  we  have 
the  energy  of  the  signal  £  =  11.284  x  10“®.  The  noise  spectral  density  Nq  = 
2.2568  X  10“®  so  that  the  energy  to  noise  ratio  is  =  10.  As  an  example 

assume,  the  gain  at  the  sensors  is  taken  as  G  =  60  dB.  For  this  setup  we  have  the 
FIM  as 

_  r  23.7315  17.8134  ' 

1-lr+l^-  ^  20.1021 

where 


■  2.9743 

2.1807  ■ 

and  X^  = 

■  20.7572 

15.6327  ■ 

2.1807 

2.3558 

15.6327 

17.7463 

(41) 


Now,  with  the  same  setup  but  with  a  signal  of  RMS  bandwidth 
the  FIMs  were 


7456.6  5467.3 
5467.3  5907.3 


5.55  MHz 
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Physical  Setup  Transmitted  signal 


Figure  21.  Physical  setup  and  Transmitted  signal 


where 


■  7435.9 

5451.7  ■ 

and  = 

■  20.7572 

15.6327  ■ 

5451.7 

5889.5 

15.6327 

17.7463 

(42) 


As  the  bandwidth  of  the  signal  increased,  the  information  in  the  time-delays 
increased  but  the  information  in  the  azimuth  modulation  remained  the  same. 
□ 

Therefore,  for  signals  with  RMS  bandwidth  in  the  order  of  100  kHz  and  ranges 
in  the  order  of  a  couple  of  kilometers,  the  information  in  the  azimuth  modulation 
is  comparable  to  the  information  in  the  time-delays.  For  larger  distance  or  for 
signals  with  larger  bandwidths,  the  information  in  the  time-delays  far  exceeds  the 
information  in  the  azimuth  modulation  and  so  even  if  there  is  any  increase  in 
azimuth  modulation  due  to  presence  of  buildings  it  may  not  signihcantly  increase 
the  total  information. 

5.3.2  With  an  Obstacle 

Next,  suppose  that  there  is  an  obstacle.  Some  of  the  receivers  may  not  have 
a  direct  line  of  sight  to  the  emitter.  To  analyze  one  aspect  of  the  obstacles,  viz. 
blocking  of  the  signal,  let  us  assume  that  they  do  not  reflect  the  signal.  So,  there 
is  no  multipath  and  so  the  signal  received  at  the  sensors  is  still  given  by  (38).  But 
now  the  AjS  are  affected  by  the  location  of  the  obstacles.  For  simplicity  assume 
that  there  is  one  obstacle  at  (x^,|/^).  This  obstacle  induces  some  kind  of  azimuth 
modulation,  say  fg{x^,y,^)  so  that  we  have 

^  ^  ^  f Vt) 
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Now,  the  information  matrix  is  given  by  (  see  Appendix  ) 


KyXj,  ,  Uj,  ) 


(Nom 


cos^  sin  'ipi  cos 

sin  xfji  cos  "04  sin^ 


(43) 

where  the  matrix  Jj  is  given  in  (D.3).  Compare  this  to  (40)  and  notice  that  the 
information  matrix  is  now  scaled  by  a  factor  of  This  factor  is  usually 

less  than  one  since  the  obstacles  do  not  induce  any  gain.  Also,  notice  that  there 

is  an  additional  term  This  term  is  almost  zero  everywhere 


except  around  the  lines  joining  the  obstacle  to  each  of  the  sensors.  This  can  be 
seen  in  Figure  22.  This  is  a  plot  of  the  azimuth  modulation  at  the  sensor  at 
(250,  0)m  (  magenta  color  )  as  a  function  of  the  emitter  location.  Notice  that 
as  the  emitter  moves  away  from  the  line  joining  the  sensor  and  the  obstacle,  the 
attenuation  decreases.  The  derivatives  of  the  azimuth  modulation  function  around 
these  lines  are  signihcant  and  so  the  contribution  of  the  Jj  matrix  to  the  total 
information  matrix  increases. 


A.(x  ,y  ) 
1  T  T 


Figure  22.  The  azimuth  modulation  function,  Ai[x^,y^)  from  (44),  at  the  sensor 
at  (250,  0)m  when  there  is  an  obstacle  at  (600,  400)m. 
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5.3.3  Example  2  -  A  Building  as  an  Obstacle 

Now  suppose  that  there  is  a  building  at  (600,400)m  as  shown  in  Figure  23. 
The  building  is  going  to  block  the  signal  if  it  is  exactly  between  the  emitter  and 


Physical  Setup 


Figure  23.  Setup  with  obstacle  very  close  to  the  line  of  sight  from  the  emitter  to 
the  sensor  at  (250,  0)m 


the  sensor.  We  model  the  blocking  of  the  signal  by  the  building  as  a  Gaussian 
pulse,  i.e  the  signal  at  a  sensor  is  completely  blocked  when  the  building  is  exactly 
between  the  sensor  and  the  emitter  and  the  attenuation  loss  reduces  as  a  Gaussian 
pulse  on  either  side  perpendicular  to  the  direct  line  of  sight.  The  attenuation  also 
reduces  as  the  distance  of  the  emitter  from  the  building  increases.  So  we  have, 


1  -  exp  ^Q)z^ 


(44) 


where  z  =  [{x^ 
given  by 


Xg)  {Vj,  —  The  matrix  C  and  the  rotation  matrix  Q  are 

cos  9  sin  9 
—  sin  9  cos  9 


C  = 


0 


a 


D 


Q  = 


where  and  determine  the  width  of  the  building  and  how  the  attenuation 
changes  with  distance  from  the  building.  The  angle  9  is  the  angle  of  the  line 
joining  the  building  and  the  sensor  i  from  the  positive  x-axis  which  is  given  by 

9  =  tan“^  I  — — ^  I .  The  rotation  matrix  is  required  to  orient  the  maior  axis 
\xi-XgJ 

of  the  Gaussian  bell  along  the  line  joining  the  building  and  the  sensor.  Figure 
22  shows  this  azimuth  modulation  at  the  sensor  at  (250,  0)m  as  a  function  of  the 
target  position  when  a  building  is  at  (600, 400)m. 

Using  the  same  setup  as  in  Example  1  but  now  with  the  building  placed  at 
(600,400)m  as  in  Figure  23  there  was  a  signihcant  increase  in  the  FIM.  The  new 


41 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED. 


FIM  was  (  compare  this  to  (41)  ) 


299.4020  -352.6595 

-352.6595  416.8441 


where 


■  0.3661 

0.2684  ■ 

and  = 

299.0359 

-352.9279  ' 

0.2684 

0.2900 

-352.9279 

416.5541 

Figure  24  shows  the  localization  ellipses  for  the  two  information  matrices.  The 
green  represents  the  localization  ellipse  when  the  building  is  not  present  and  the 
red  represents  the  localization  ellipse  when  the  building  is  present.  Notice  that 
the  minor  axis  is  considerably  reduced  while  the  major  axis  is  elongated.  This  is 
because  when  the  building  blocks  the  signal  to  a  sensor,  we  know  that  the  emitter 
is  on  the  line  joining  the  building  and  the  sensor.  Here  the  gain  G  was  reduced  to 
5000  to  make  the  ellipses  £t  in  the  figure. 


Comparision  of  Information  using  localization  ellipses 


Figure  24.  Localization  ellipses  for  the  information  matrices.  Green  -  no  building. 
Red  -  With  building. 


Now,  when  the  building  is  placed  at  the  point  (420,  400)m  which  is  away  from 
the  lines  joining  the  emitter  to  the  sensors  as  shown  in  Figure  25,  the  information 
matrix  was  not  affected  and  remained  the  same  as  in  Example  1. 

Another  interesting  fact  is  that,  for  the  azimuth  modulation  function  we  have 
assumed  in  (44),  the  derivative  at  any  point  exactly  on  the  line  joining  the  sensor 
and  building  is  less  than  at  any  point  on  either  side  of  the  line.  So,  when  we  place 
the  building  exactly  on  that  line  at  (600,  407.9137)m,  the  information  matrix 
reduced  to 


7.2817  9.0970 

9.0970  12.4407 
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Physical  Setup 


Figure  25.  Setup  with  the  obstacle  away  from  the  lines  of  sight  from  the  emitter 
to  the  sensors. 


Next,  with  the  same  setup  but  with  a  signal  of  RMS  bandwidth  =  5.55  MHz 
the  FIMs  were  (  compare  this  to  (42)  ) 


1214.4  318.2 

318.2  1141.5 


where 


■  915.3456 

671.0949  ■ 

and  = 

299.0359 

-352.9279  ' 

671.0949 

724.9918 

-352.9279 

416.5541 

Figure  26  shows  the  corresponding  localization  ellipses.  So,  the  total  information 
is  reduced.  Even  though  the  building  increased  the  information  in  the  azimuth 
modulation  it  reduced  the  information  in  the  time-delays  resulting  in  a  net  decrease 
in  the  total  information.  □ 


5.4  Conclusions 

Under  certain  specihc  circumstances  having  the  building  can  increase  the 
Fisher  information  matrix  if  the  building  is  close  to  the  lines  of  sight  from  the 
emitter  to  the  sensors.  This  is  because,  when  the  signal  at  any  one  sensor  is 
considerably  attenuated  then  we  know  that  the  emitter  is  on  the  line  joining  the 
building  and  the  sensor.  This  narrows  down  the  possible  emitter  position  from  a 
plane  to  a  single  line.  Then  the  signals  received  at  the  other  sensors  are  useful  in 
determining  the  precise  location  of  the  emitter  on  that  line.  When  the  building  is 
not  close  to  the  direct  path  between  the  sensor  and  the  emitter  then  there  is  no 
increase  in  information  due  to  the  presence  of  the  building. 
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Physical  Setup 


Figure  26.  Localization  ellipses  for  the  information  matrices.  Green  -  no  building. 
Red  -  With  building. 
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6  Improved  TDOA  Position  Fixing 
6.1  Introduction 

Localization  is  the  process  of  extracting  the  location  information  of  a  signal 
transmitter  by  intercepting  the  transmitted  signal.  The  location  information  can 
be  obtained  from  one  or  more  of  the  received  signal  properties  such  as  signal 
strength,  angle  of  arrival,  TDOA  etc.  In  this  paper  we  focus  on  the  TDOA  based 
localization,  also  called  multilateration.  A  conventional  multilateration  has  two 
steps.  In  the  first  step  the  TDOAs  between  at  least  four  pairs  of  sensors  are 
estimated.  In  the  second  step,  these  TDOAs  are  used  along  with  the  known  sensor 
locations  for  position  hxing.  Abundant  literature  exists  focusing  on  one  or  both 
of  these  steps.  Much  of  the  original  work,  particularly  in  sonar,  modeled  the 
transmitted  signal  as  a  stochastic  process  with  the  power  spectrum  known,  in 
active  localization  or  unknown,  in  passive  localization  [25,  42,  24,  21].  The  noise 
at  each  sensor  is  almost  always  assumed  to  have  a  known  power  spectrum  and 
uncorrelated  from  sensor  to  sensor.  For  this  model  and  assumptions,  Hahn  and 
Tretter  [42]  proposed  a  delay  vector  estimator  using  correlators.  They  showed 
that  when  there  are  M(>  2)  sensors,  even  though  we  are  only  interested  in  the 
M  —  1  TDOAs  with  respect  to  a  reference  sensor,  it  is  benehcial  to  compute 
the  M(M  —  l)/2  TDOAs  for  all  possible  sensor  pairs  combinations  using  cross¬ 
correlation  and  then  use  the  Gauss-Markov  estimate  of  the  desired  M  —  1  TDOAs 
with  respect  to  the  reference  sensor.  They  also  showed  that  such  TDOA  estimates 
attain  the  CRLB.  Knapp  and  Carter  [21]  on  the  other  hand  proposed  a  generalized 
correlation  method  for  the  estimation  of  the  TDOAs  where  in,  they  pre-hlter  the 
received  signals  in  order  to  maximize  the  SNR  before  correlating.  Fowler  and  Hu 
have  shown  that  the  above  model  and  assumptions  are  reasonable  for  sonar  but 
not  for  radar  [28].  They  have  showed  that,  for  passive  localization  using  radars, 
the  signal  must  be  modeled  as  deterministic  and  unknown.  For  this  signal  model 
assumption  Stein  [9]  has  derived  the  MLE  for  the  TDOA  between  two  sensors. 
Stein’s  result,  similar  to  Hahn  and  Tretter’s  [42],  is  that  the  MLE  for  the  TDOA 
between  two  sensors  is  the  differential  delay  that  maximizes  the  cross-correlation 
of  the  signals  received  at  the  two  sensors. 

Irrespective  of  the  signal  model  assumptions  and  how  the  TDOAs  are  esti¬ 
mated,  each  TDOA  between  a  pair  of  sensors  “fixes”  the  position  of  the  emitter  on 
a  hyperboloid  with  foci  at  the  locations  of  the  two  sensors.  In  2-dimensions,  a  set  of 
three  TDOAs  uniquely  dehne  the  location  of  the  emitter  as  the  single  intersection 
point  of  their  respective  hyperbolas.  Therefore,  position  hxing  using  the  TDOAs  is 
a  non-linear  problem.  Torrieri  [17]  proposed  a  linear  least  squares  estimator  where 
this  nonlinear  relation  between  the  TDOAs  and  the  emitter  location  is  linearized 
by  expanding  it  using  Taylor  series  about  a  reference  point  and  retaining  the  hrst 
two  terms.  Ho  and  Chan  [43,  18]  derived  a  closed  form  two-step  WLS  estimator 
that  is  asymptotically  efficient.  They  hrst  introduce  an  intermediate  variable  that 
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is  a  function  of  the  emitter  location  to  linearize  the  nonlinear  equations  and  hnd 
the  least  squares  (LS)  estimator.  Next,  they  use  the  relation  between  the  emitter 
location  and  the  intermediate  variable  to  solve  a  second  WLS  to  arrive  at  the  hnal 
solution.  Recently,  Matthew  and  Silverman  [20]  have  derived  a  simple  closed-form 
least  square  estimator  that  performs  quite  well. 

Alternately,  when  the  signal  waveform  is  known,  as  in  active  localization  for 
example,  the  received  signals  are  correlated  with  the  known  signal  waveform  to 
estimate  the  TOAs  [39].  Here,  it  is  common  to  work  with  the  TOAs  instead  of 
the  TDOAs  [29,  30,  44].  The  TOA  at  a  sensor  “hxes”  the  emitter  position  on 
a  sphere  (circle  in  two  dimensions)  centered  at  the  sensor  location.  Here,  the 
intersection  point  of  the  spheres  is  taken  as  the  estimate  of  the  emitter  location. 
In  the  TOA  based  estimation  the  unknown  transmission  time  often  occurs  as  a 
nuisance  parameter  that  needs  to  be  estimated  as  well.  Although,  Do  et.  al 
[31]  have  shown  that  the  TOA  and  TDOA  measurements  are  transformable  to 
each  other  without  loss  of  information  regarding  positioning  and  thns  the  position 
estimators  based  on  them  should  be  theoretically  equivalent. 

These  position  hxing  techniqnes  are  good,  even  efficient,  in  nsing  the  infor¬ 
mation  in  the  TDOAs.  Bnt  there  is  more  information  in  the  received  signals  than 
jnst  the  TDOAs.  For  example,  a  measnre  of  the  accnracy  of  each  of  the  TDOAs 
can  also  be  obtained.  This  measnre  can  be  used  to  weight  the  TDOAs  so  that  the 
more  accurate  TDOA  estimates  have  a  larger  contribution  in  the  position  fixing 
step.  Snch  WLS  type  techniqnes  have  also  been  proposed  in  the  literatnre.  Bnt 
often  the  weighting  matrix  is  the  inverse  of  the  covariance  matrix  of  the  TDOAs 
which  is  conveniently  assnmed  to  be  known  [27,  19].  In  some  cases  where  the  SNR 
is  known,  it  is  used  for  the  weighting. 

Recently,  Yeredor  and  Angel  [23]  have  proposed  a  method  where  the  FIM  is 
used  as  the  weighting  matrix.  In  order  to  compnte  the  FIM,  the  emitter  location 
must  be  known.  So  they  nse  an  iterative  techniqne  where  they  perform  estimation 
over  consecntive  observation  intervals  and  nse  the  emitter  location  estimated  with 
the  previons  observation  interval  for  compnting  the  FIM  to  be  used  in  the  current 
observation  interval.  Also,  the  CRLB  (  the  inverse  the  FIM  )  is  “on  the  average” 
the  asymptotic  variance  of  the  MLE.  So,  it  is  not  a  very  accnrate  measnre  of  the 
qnality  of  the  TDOA  estimate  from  that  particular  observation.  In  this  paper  we 
propose  a  novel  approach  to  estimating  the  covariance  matrix  simnltaneously  with 
the  TDOAs  from  a  given  observation  interval  and  then  use  it  for  the  position  hxing 
step  to  improve  the  localization  performance. 

6.2  Methods,  Assumptions,  and  Procedures 

Onr  focus  is  on  radar  and  so  we  use  the  signal  modeling  assnmptions  of 
Stein  [9].  Snppose  that  a  stationary  emitter  is  located  at  an  nnknown  location 
and  a  network  of  M  sensors  are  located  at  known  locations  i  = 

0, 1,  •  •  •  ,  M  —  1  as  shown  in  Fignre  27.  For  simplicity  we  are  assnming  a  two  di¬ 
mensional  case.  Extension  to  the  three  dimensional  case  is  straightforward.  The 
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Physical  Setup 
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Figure  27.  Physical  placement  of  the  sensors  (for  M=4)  and  the  Emitter  position 
used  for  simulation. 


sensors  are  all  synchronized  in  time  and  each  of  the  sensors  intercepts  the  signal 
within  the  time  interval  (0,T).  The  emitter  transmits  an  unknown  signal  s{t)  for 
an  unknown  duration  Tg  <  T  starting  at  an  unknown  time  to  <  T.  We  shall  as¬ 
sume  that  the  transmitted  signal  s{t)  is  real.  It  can  be  narrowband  or  wideband, 
lowpass  or  bandpass.  After  interception,  the  signal  received  at  sensor  i  in  the 
presence  of  noise  can  be  written  as 


ri(t)  =  Ais(t  -  Ti)  -f  Wiit),  0  <  t  <  T,  ,  . 

where  Wi{t)  is  a  zero  mean  wide  sense  stationary  additive  white  Gaussian  random 
process,  Aj’s  are  the  unknown  attenuations  due  to  propagation  loss,  assumed  real, 
and  the  r^’s  are  the  unknown  TOAs  given  by 


Ti  = 


-XiY  -h  (i/j,  -ViY 


+  to,  i  =  0,  -  ■  ■  ,M  -  1 


(46) 


where  c  is  the  propagation  speed  of  the  signal.  We  assume  that  the  noise  at  a 
sensor  is  independent  of  the  noise  at  any  other  sensor,  i.e.,  Wi(t)  and  Wj(t)  are 
independent  for  i  j. 


6.3  Variance  of  the  TDOAs 

Let  the  vector  of  the  TOAs  be 

T  =  [tq  •  •  • 

and  the  vector  of  TDOAs  with  respect  to  sensor  0  be 


(47) 


(48) 
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where  t-  =  Ti  —  Tq,  i  =  1,  ■■■,  M  —  1  and  (=t:)^  represents  the  transpose  of  the 
matrix  (*).  Stein  has  derived  the  MLE  for  the  TDOA  between  a  pair  of  sensors. 
The  MLE  for  r/  is  given  by 

f-  =  argmax  /  ri(t)ro(r +  t)  dt.  (49) 

Jo 


In  order  to  estimate  the  TDOAs  at  all  sensors  with  respect  to  the  reference  sensor  a 
common  practice  is  to  use  Stein’s  MLE  to  estimate  each  of  the  TDOAs  individually. 
That  is, 

•••  (50) 

Note  that  this  (M  —  1)  x  1  TDOA  vector  estimate  in  (50)  is  not  the  MLE  of 
the  TDOA  vector  in  (48).  In  fact,  when  M  >  2,  Hahn  and  Tretter  [42]  showed 
that  computing  the  M(M  —  l)/2  TDOAs  for  all  possible  sensor  pairs  using  cross¬ 
correlation  and  then  using  the  Gauss-Markov  estimate  of  the  desired  M  —  1  x  1 
vector  of  TDOAs  with  respect  to  the  reference  sensor  is  better  than  (50)  and 
attains  the  CRLB.  When  M  is  large  computing  all  the  M(M  —  l)/2  TDOAs  may 
be  impractical.  We  will  continue  to  use  (50)  as  the  estimator  for  (48)  as  our  focus 
is  on  improving  the  localization  accuracy  by  weighting  the  TDOAs.  Instead,  we 
will  further  use  (49)  to  also  estimate  the  variances  of  each  of  the  f-s. 

It  is  a  well  known  fact  that  the  argument  that  maximizes  the  likelihood  func¬ 
tion  is  the  MLE  while  the  asymptotic  variance  is  equal  to  the  negative  of  the 
expected  value  of  the  curvature  of  the  likelihood  function  at  the  peak.  Therefore, 
we  will  use  the  eurvature  of  the  likelihood  function  as  the  estimate  of  the  inverse  of 
the  variance  of  the  corresponding  TDOA.  This  indicates  the  quality  of  that  particu¬ 
lar  estimate  and  not  the  “on  the  average”  type  quality  measure  as  with  the  CRLB. 
The  likelihood  function  in  (49)  is  the  correlation  function  itself  given  by 


R{t)  =  [  r-iitypij +  t)  dt.  (51) 

Jo 

Figure  28  shows  segments  of  two  correlation  functions  around  their  peak,  for  dif¬ 
ferent  SNRs.  Here  the  SNR  at  the  reference  sensor  is  set  at  0  dB.  The  SNRs  at 
two  other  sensors  that  have  the  same  TDOA  are  set  at  0  dB  and  —10  dB.  The 
signals  received  at  the  two  sensors  are  correlated  with  the  signal  at  the  reference 
sensor.  The  peak  value  is  subtracted  from  the  corresponding  correlation  function 
to  bring  them  to  the  same  level  for  visual  comparison.  We  dehne  the  SNR  at 
each  sensor  as  the  average  power  of  the  received  signal  to  the  average  noise  power 
at  the  sensor,  i.e,  if  at  the  ith  sensor,  Vsi  =  /g  |’S(t)p  dt  is  the  average 

power  of  the  signal  and  ^  is  the  noise  spectral  density,  then  the  SNR  is  given  by 

10  log  ^ /2^^x  BW^  where  BW  is  the  bandwidth  of  the  receiver.  Notice 

that  the  signal  with  the  higher  SNR  has  a  larger  curvature  indicating  a  more  accu¬ 
rate  TDOA  and  hence  should  be  weighted  more  heavily.  Now,  the  actual  curvature 
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Figure  28.  Curvatures  of  the  correlation  function  at  the  peak  for  different  SNRs. 

of  the  likelihood  function  depends  on  both  the  transmitted  signal  waveform  s{t) 
and  the  SNR,  both  of  which  are  unknown.  So,  we  need  to  use  some  curve  htting 
techniques  around  the  peak  to  estimate  the  curvature  of  the  likelihood  function. 
We  use  the  simplest  form  of  htting  viz.  quadratic  htting  around  the  peak. 


6.4  Results  and  Discussion 

For  the  simulation,  we  used  the  setup  shown  in  Figure  27.  We  used  a  Gaussian 
chirp  as  the  unknown  transmitted  signal  which  is  given  by 


s(t)  =  exp 


sin(27rmt^) 


(52) 


The  signal  is  assumed  to  be  approximately  zero  for  t  <  0  and  for  t  >  Tg. 
We  set  Tg  =  5/is  and  =  27i  x  10^.  The  observation  interval  at  each  of  the 
sensors  was  taken  to  be  T  =  0.2  ms.  The  unknown  transmission  time  of  the 
signal  was  set  to  to  =  0.07  ms.  With  this  conhguration  the  maximum  TDOA  is 
0.0988  ms.  The  frequency  spectrum  of  the  Gaussian  window  is  given  by  |5'(F)|  = 
(\/^/(7j,,)  exp  (— 27r^F^/(j^)  [39],  and  the  bandwidth  of  the  chirp  is  Bg  =  1.5 
MHz.  The  rate  of  change  of  frequency  for  the  linear  chirp  was  chosen  to  be  m  = 

7^  =  3  X  10®.  The  correlation  integrals  are  approximated  using  summations  with 

S 

5t  =  0.33  ns  which  means  the  sampling  frequency  is  Fg  =  300  MHz.  The  reason 
for  choosing  such  high  sampling  frequency  is  that  at  this  frequency,  the  position 
quantization  error  due  to  sampling  is  in  the  order  of  {c/ Fg)  =  10“®  km. 

The  regular  two-step  LS  TDOA  approach  is  implemented  as  follows.  In  the 
hrst  step,  the  TDOA  f/  estimates  were  obtained  by  cross-correlating  the  signal  at 
each  sensor  with  the  signal  at  sensor  0.  Then  a  1km  x  1km  region  around  the  true 
emitter  location  was  split  into  10  x  10  grid  points  and  for  each  emitter  location  on 
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the  grid  point  the  TDOAs  were  computed  using  the  formula 


for  i  =  1,2,  ••  •  ,M  —  1.  Next  the  LSE  between  the  estimated  TDOAs  and  the 
computed  TDOAs  was  calculated  as 

M-l 

LSE  =  5;  (f'  -  r')^  (54) 

i=l 

This  LSE  is  a  function  of  the  emitter  location  The  emitter  location  that 

minimized  the  LSE  is  the  LS  estimate  of  the  emitter  location.  Since  we  know 
that  the  LSE  is  a  2-dimensional  parabolic  function  of  the  emitter  location  near  the 
minimum,  we  improved  the  accuracy  by  htting  a  parabola  through  the  100  points 
at  which  the  LSE  was  computed.  The  WLS  estimator  is  also  implemented  in  the 
same  manner  except  the  WLS  error  (WLSE)  was  calculated  as 

M-l 

WLSE  =  ^  aM  -  T[f  (55) 

i=l 

where  af’s  are  the  weights.  The  emitter  location  that  minimized  the  WLSE  is  the 
WLS  estimate  of  the  emitter  location.  In  order  to  compute  the  weights,  we  chose  a 
window  of  30  points  around  the  true  peak  of  the  correlation  function.  Using  these 
30  points  we  computed  the  coefficients  (pi,P2,P3)  of  a  parabola  +  p2t  +  ps 
that  hts  these  points  in  the  least  squares  sense.  Then  the  absolute  value  of  the 
curvature  of  the  parabola  was  taken  as  the  weights,  i.e, 

a=|2pi|  (56) 

We  ran  1000  Monte  Carlo  simulations  to  compare  the  performance  of  the  two 
estimators.  The  SNR  was  set  at  0  dB  at  the  reference  sensor  and  at  the  other  two 
sensors.  At  the  fourth  sensor,  the  SNR  was  varied  from  —30  dB  to  —20  dB.  Figure 
29  shows  the  comparison  of  the  mean  square  error  of  the  1000  estimates  for  each 
of  the  two  estimators  for  different  SNR  values  of  the  fourth  sensor. 

Notice  that  the  performance  of  the  LS  TDOA  approach  deteriorates  signif¬ 
icantly  as  the  SNR  goes  below  —28  dB.  This  is  because  at  such  low  SNRs  the 
TDOA  estimator  in  (49)  completely  breaks  down.  When  the  noise  level  is  low,  the 
peak  location  of  (51)  is  close  to  the  true  TDOA  value  but  when  the  noise  level  is 
sufficiently  high,  the  peak  of  (51)  can  occur  anywhere  in  (— T,  T)  with  high  prob¬ 
ability.  Figure  30  shows  the  histogram  of  the  peak  of  the  correlation  function  in 
(51)  for  the  SNR  values  of  0  dB  and  —30  dB.  Notice  that  at  0  dB  the  peak  location 
lies  very  close  to  the  true  TDOA  value  of  —0.0987  ms  while  at  —30  dB,  the  peak 
of  a  realization  appears  at  —0.06  ms  also.  When  this  outlier  is  used  in  the  second 
step  for  position  hxing,  the  location  estimate  is  far  from  the  true  location.  In  the 
WLS  TDOA  approach,  the  effect  of  this  outlier  TDOA  estimate  is  minimized  by 
the  weighting. 
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Figure  29.  MSE  of  the  1000  estimates  for  each  of  the  two  estimators  for  different 
SNR  values  of  the  fourth  sensor. 
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Figure  30.  Histograms  of  the  peak  location:  (a)  SNR  =  0  dB.  (b)  SNR  =  —30  dB. 

6.5  Conclusion 

We  have  proposed  a  weighted  least  squares  position  fixing  technique  for  the 
multilateration  problem.  We  showed  that  the  TDOA  estimator  breaks  down  for 
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very  low  signal-to-noise  ratios  and  when  such  TDOA  estimates  are  used  for  position 
hxing,  the  location  estimates  are  far  from  the  true  emitter  location.  Our  weight¬ 
ing  approach  mitigates  the  effect  of  such  outliers  in  the  position  hxing  step.  We 
proposed  a  simple  technique  for  the  computation  of  the  weights  from  the  received 
signal. 
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7  Overall  Conclusion 


We  have  addressed  the  problem  of  passively  gathering  information  about  an 
emitter  of  electronic  signals  using  distributed  sensors.  First  we  proposed  an  asymp¬ 
totically  optimal  technique  for  detection  of  the  presence  or  absence  of  such  signals 
in  the  data  collected  at  the  distributed  sensors.  This  is  a  centralized  detector 
unlike  the  commonly  used  decentralized  decision  fusion  techniques.  The  derived 
detector  is  the  generalized  likelihood  ratio  test  (GLRT)  detector.  We  also  derived 
simpler  detectors  from  the  GLRT  by  making  various  assumptions.  Receiver  oper¬ 
ating  characteristics  (ROG)  curves  for  currently  used  detectors  are  computed  and 
compared  to  the  ROG  curves  for  the  GLRT  detector.  After  the  presence  of  such 
signals  is  detected,  the  next  step  is  to  estimate  the  location  of  the  emitter.  For 
this  we  proposed  the  maximum  likelihood  estimator  (MLE).  The  conventional  ap¬ 
proach  for  localization  using  multiple  sensors  is  to  hrst  estimate  the  time  difference 
of  arrivals  (TDOAs)  of  the  signals,  independently  between  pairs  of  sensors  and  then 
to  hnd  the  location  of  the  emitter  using  the  intersection  point  of  the  hyperbolas 
dehned  by  these  TDOAs.  This  is  referred  to  as  the  conventional  TDOA  technique 
and  it  has  been  shown  in  the  literature  that  this  two-step  approach  is  suboptimal 
in  comparison  to  what  is  called  the  direct  position  determination  (DPD)  approach. 
In  the  DPD  approach,  the  intermediate  step  of  estimating  the  TDOAs  is  bypassed 
and  the  location  is  estimated  directly  from  the  observations.  In  this  paper  we 
take  the  DPD  approach  instead  of  the  conventional  two-step  approach.  The  DPD 
type  localizers  that  have  been  proposed  in  the  literature  are  based  on  certain  as¬ 
sumptions  on  the  transmitted  signal  such  as  narrowband  or  wideband,  lowpass  or 
bandpass  etc.  We  make  no  such  assumptions  on  the  signal  and  this  paper  covers 
a  wide  variety  of  transmitted  signals.  In  passive  localization,  it  is  common  to  not 
know  the  transmission  time  of  the  signal  and  more  often  the  signal  waveform  it¬ 
self  is  unknown.  So,  we  have  analyzed  these  two  commonly  occurring  cases  of  (i) 
signal  waveform  unknown  and  (ii)  signal  waveform  known  with  unknown  trans¬ 
mission  time.  The  localizers  proposed  in  literature  assumed  discrete  time  but  they 
have  not  addressed  the  quantization  like  effect  on  the  location  estimate  due  to 
sampling  of  the  received  signals.  To  avoid  this  quantization  like  effect,  we  have 
used  a  continuous  time  model.  We  have  also  derived  the  Fisher  Information  Matrix 
(FIM)  which  gives  a  deeper  insight  into  the  relations  between  various  parameters 
in  the  model  and  their  identihability.  We  showed  that  the  proposed  MLE  outper¬ 
forms  the  conventional  two-step  localizers  and  also  attains  the  Gramer  Rao  Lower 
Bound  (GRLB)  for  high  signal-to- noise  ratios  (SNR).  Though  the  performance  of 
the  MLE  attains  the  GRLB,  there  is  still  scope  for  improvement.  This  improve¬ 
ment  comes  from  the  geometry  of  the  sensors.  We  showed  that  for  a  given  SNR, 
the  GRLB  depends  on  the  sensor  geometry.  Thus,  the  sensor  geometry  can  be 
optimized  in  order  to  further  reduce  the  GRLB  and  also  the  variance  of  the  MLE. 
We  have  dehned  this  problem  of  optimizing  the  sensor  geometry  and  derived  the 
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optimal  sensor  configurations  for  a  few  specific  scenarios.  The  optimization  of  the 
sensor  geometry  for  a  general  scenario  is  quite  elusive  and  still  open  for  research. 
In  addition  to  efficient  signal  processing  as  in  the  MLE  and  optimizing  the  sensor 
geometries,  the  localization  performance  can  further  be  enhanced  by  using  infor¬ 
mation  about  the  terrain.  This  is  commonly  referred  to  as  the  knowledge  aided 
design  (KAD).  KAD  is  particularly  useful  for  localization  in  urban  environments 
where  the  distances  are  small  and  the  knowledge  of  the  terrain  is  very  well  known. 
We  investigated  the  aspect  of  azimuth  modulation  induced  by  various  objects  in 
the  terrain  and  the  usefulness  of  this  modulation  for  localization.  In  particular, 
we  have  shown  that  when  an  obstacle  blocks  the  signal  to  one  of  the  sensors,  then 
there  is  increase  in  the  over  all  information  of  the  location  of  the  emitter.  All  the 
concepts  mentioned  above  are  for  improving  the  localization  performance  and  do 
require  altering/increasing  the  existing  physical  resources.  For  example  the  MLE 
requires  high  bandwidth  links  between  all  the  sensors  and  the  fusion  center  in  order 
to  transmit  all  the  data  collected  at  each  sensor  for  simultaneous  processing  at  the 
fusion  center.  This  is  in  contrast  to  the  conventional  TDOA  which  only  requires 
low  bandwidth  links  to  the  fusion  center  because,  here  only  the  TDOA  instead  of 
the  complete  observation,  is  transmitted  to  the  fusion  center.  So  hnally  we  pro¬ 
pose  an  improvement  to  the  conventional  TDOA  approach  which  does  not  require 
any  additional  physical  resources  but  still  signihcantly  increases  the  localization 
performance  particularly  for  SNRs  at  the  break  down  range.  The  only  additional 
piece  of  information  that  needs  to  be  transmitted  to  the  fusion  center  along  with 
the  TDOA  is  the  curvature  of  the  likelihood  function.  In  the  conventional  TDOA 
approach  the  hrst  step  is  to  estimate  the  TDOAs  and  the  second  step,  called  the 
position  hxing,  is  to  estimate  the  emitter  location  as  the  intersection  of  the  hyper¬ 
boloids  dehned  by  these  TDOAs.  For  the  TDOA  estimation  the  commonly  used 
estimator  is  the  time-delay  that  maximizes  the  cross-correlation  function.  This 
is  the  MLE  of  the  TDOA  (  note  that  this  is  not  the  MLE  of  the  emitter  loca¬ 
tion  which  we  have  derived  )  and  the  cross-correlation  function  is  the  maximum 
likelihood  function.  Now,  since  the  asymptotic  variance  of  an  MLE  is  equal  to 
the  negative  of  the  expected  value  of  the  curvature  of  the  likelihood  function,  we 
proposed  a  weighted  least  squares  type  position  hxing  technique  where  the  weights 
can  be  computed  from  the  curvature  of  the  likelihood  function.  Hence,  we  have 
addressed  the  problem  of  passively  gathering  information  about  an  emitter  of  elec¬ 
tronic  signals  using  distributed  sensors  and  investigated  various  aspects  that  can 
increase  this  information. 
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A  Derivation  of  the  GLRT  Detector 


The  GLRT  decides  "Hi  if 


p(r;s,  A,h,k,?^i) 

or  equivalently  if  the  test  statistic 

InLG(f)  >  Iny  =  7' 

Here  s,  A,h  and  k  are  the  MLEs  of  s,  A,  n  and  k  respectively  assnming  "Hi  is 
true.  Since  w  ~  CA/'(0,  Imv)  we  have 


LrAr)  = 


^exp  (  -  (r  -  H(A,n,k)s)  (r-H(A,n,k)s 


^  exp  (-f^f) 


Taking  logarithms  we  have 


InLG'(r)  =  —  (  r  —  H(A,  h,  k)s )  (  r  —  H(A,  h,  k)s  )  +  r^f 


2Re  f'^H(A,  h,  k)s  —  ^H(A,h,  k)sj  ^H(A,h,  k)s 


The  MLEs  s,  A,  h  and  k  maximize  p(f;  s,  A,  n,  k;  "Hi).  But, 

p(f;s,  A,n,k,?^i)  =  — ^exp  -  ff  -  H(A,  n,  k)s)  ff-H(A,n,k)s 


Maximizing  p(r;  s,  A,  n,  k,  "Hi)  over  s,  A,  n  and  k  is  equivalent  to  minimizing  the 
exponent  which  is 


J  =  ( f  —  H(A,  n,  k)s )  (r  — H(A,n,  k)s 


From  [39]  we  have  the  MLE  of  s  as 

s=  H^(A,  n,  k)H(A,  n,  k)  H(A,  n,  k)f 
Putting  this  back  in  the  test  statistic  we  have 

InLG(f)  =  max  f^H(A,  n,  k)  [h^(A,  n,  k)H(A,  n,  k)!  ^  H^(A,  n,  k)i 

A,n,k  I  L 
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Notice  that 


H  (A,n,k)H(A,n,k)  =  (A  A)l7v 


H^(A,n,k)f  =  (  ^ 


so  that  we  have 


In  LMi)  =  max 


r^H(A,n,  k)H  (A,n,  k)f 


=  max 

A,n,k 


H^(A,n,  k)rj  ^H^(A,  n,  k)r 
(A^A) 


=  max 

A,n,k 


(A^A) 


=  max 

A,n,k 


\  /M-1 

jk,  ^ 

A*  (■Dni\H  z 


Y,  Aff  (P"0(W  ‘)  A*(W 


(A^A) 


M-l  M-1 


=  max 

A,n,k 


i=0  j=0 


(A^A) 


If  we  let  B(n,  k)  be  the  cross  ambiguity  matrix  dehned  as 


B(n,k)  = 


f^pnow"“ 


fyp"ow'“0(w'“0)®(P"0)-H‘i;Q 

fyp"iw'‘i(w'“o)®(p”o)-fffQ 


fyp”ow*’o  (w'‘“-i)®(p"'“-i)®?M-i 
ff  p"l  W*’!  (w'‘“-l)®(P”J'^-l)®fM-l 


f  M  - 1 P  "  “  “  ^  W  - 1  ( W  0  ) '^  (  P  "■  0  ) -ff  f  0 


then  we  have  a  Hermit  ian  form  and  so 


,  ^  A^BVn,k)A  ,  ,  ,, 

In  Lcir)  =  max - - =  max  Amax(B  (n,  kj) 

A,n,k  A  A 
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where  A^ax  is  the  maximum  eigenvalue  of  B  (n,  k).  The  MLE  of  A  is  the  eigen¬ 
vector  corresponding  to  Amax(B  (n,  k)).  Therefore,  the  GLRT  decides  "Hi  if 

maxAmax(B*(n,  k))  >  7' 
n.k 

The  cross  ambiguity  matrix  B(n,  k)  is  Hermitian  with  real  and  positive  eigenvalues. 
So,  we  have  Amax(B  (n,  k))  =  Amax(B(n,  k))  and  so  the  GLRT  decides  "Hi  if 

maxAmax(B(n,  k))  >  7' 
n.k 
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B  Localization 


B.l  CRLB 

We  will  derive  the  CRLB  for  the  emitter  location  estimate.  First,  we  will  show 
that  the  FIM  for  the  model  used  for  unknown  signal  with  unknown  transmission 
time  case  is  singular.  We  will  then  use  a  transformation  of  the  parameters  in  the 
model  and  derive  the  CRLB.  Let  r  =  [tq  ri  •  ■  •  i  A  =  [Rq  Ai  •  •  • 

and  0  =  [oo  Oi  ■  ■  •  otv-i  &2  ■  ■  ■  where 

do  =  X  fo  ~  T  fo  cos27rnFot  dt,  ^  /o  <s(f)  sin  27rnFo^  dt 

Let  6  =  [r^  0^]^-  The  TOAs  r*  are  a  function  of  the  emitter  location  {x^,y^) 

and  the  signal  transmission  time  to- 


where  c  is  the  propagation  speed  of  the  signal.  If  l{6)  is  the  log-likelihood  function, 
then  the  FIM  is  given  by 


f  ] 

\  dTdT'^  j 

I  an(e) } 

[dAdT'^S 

f  an(e)  ] 

\  (90(9t^  J 


r  dHjo)  ] 

1  drdA^  J 

r  dHjO)  \ 

1  dAdA^  J 

[  dHjO)  ] 

\  (90(9 A^  J 


r  dHjo)  ]  1 
\  Ordcf)^  j 

\  (9A(90^  J 
r  9^1 

1  (90(90^  J 


(B.l) 


B.1.1  Signal  unknown  with  unknown  transmission  time 

From  (23),  we  have  the  log-likelihood  function  as 


M—l 


l{6)  =  [  Y]  {xm{t)  -  dt  (B.2) 


where  h{t)  is  as  dehned  in  Appendix  B.3.  Partial  differentiation  with  respect  to 
(w.r.t)  Ti  gives. 


L-  ~  No  I 

and  w.r.t  Aj  gives. 


dh^{t  -  Ti) 


0  ]  dt  (B-3) 


[  2  {ri{t)  -  Aih^{t  -  Ti)(t))  {-h^{t  -  Ti)(t>)  dt  (B.4) 

O/ii  iVn  ./n 
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for  i  =  0, 1,  •  •  ■  ,  M  —  1.  Partial  differentiation  w.r.t  (p  gives, 


dl{0) 

dcp 


M—l 


N, 


-^0  m=0 


2  {Xm{t)  -  Amh^{t  -  Tm)(t>)  {-Amh^{t  -  T^))  dt  (B.5) 


Next  we  evaluate  the  second  derivatives.  Partial  differentiation  of  (B.3)  w.r.t  r* 
gives, 

““il  f  (n(t)  -  A,h^(t  -  Ti)<l>)  dt 

+/J' 

Taking  the  negative  of  the  expected  value  of  both  sides  and  using  (B.27)  gives. 


f  d‘^l{0)  \  2  /  dh^it  -  n 


I  dri  /  iV, 


=  ^  /  (A, 

0  Jo 


(9r,' 


4>  )  dt 


A? 


(^^0/2) 


■  dh{t  -  T 

i)  dh^(t  -  Ti)' 

Jo  dTi 

dTi 

(p 


{T/2){2nFo)\^TrrT 


iNo/2) 


{cP^  LV  cP)  Al 


dH{G) 


where  the  2iV  — 1  x  2iV  — 1  matrix  L  is  as  dehned  in  Appendix  B.3.  Since  tt — =  0 

OTiOTj 


for  i  jtz  we  have 


-E 


drdT^ J 


(A'o/2) 


-  LL^  (/)  ( diag  ( A ) )  ^ 


(B.6) 


where  diag(A)  is  an  M  x  M  diagonal  matrix  with  Ah  diagonal  element  as  Aj. 
Partial  differentiation  of  (B.4)  w.r.t  r*  gives. 


dH{e) 

dTjdAj 


Nn 


0  L^o 


{ri(t)  -  Aih^{t  -  Ti)(p) 


dh^{t  -  Ti) 


On 


(p  I  dt+ 


^  .  dh^ it  —  Ti)  ,  \  ,  jj,  . 

-Ai - - cp  I  —  Ti)4>)  dt 


10  \  dn 

Taking  the  negative  of  the  expected  value  of  both  sides  and  using  (B.28)  gives, 

dh{t  —  Ti 


p  r  d%o)  \  A  ,T 

\  dndA,  I  (No/2) 


'*2^  I.T 


Uo 


dTi 


(t  —  Ti)  dt 


(p 


(Ao/2) 
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Since  tt — ^  =  0  for  i  ^  we  have 

OTiOAj 


dH{0)  ]  (T/2)(27rFo),^^ 


drOA^ 


{No/2) 


{(f)  L0)(diag(A)) 


(B.7) 


Partial  differentiation  of  (B.5)  w.r.t  r*  gives, 

d‘^l{6)  2  \  .  A  j  T/  \  i\  f  A  dhF {t  —  Ti)\  , 

[  ^ — —</]  {-Ah^{i  -  n))  dt 


Taking  the  negative  of  the  expected  value  of  both  sides  and  using  (B.28)  gives, 


dH{0)  ]  Af  ^ 
drM^i  (iVo/2)^ 


dh{t  -  Ti)  j, 


{t  —  Ti)  dt 


{T/2){27rFo)Af 

(iVo/2)  ^  ^ 


So,  we  have 


(T/2)(2irF„),,  „ 

(iVo/2) 


(B.8) 


Partial  differentiation  of  (B.4)  w.r.t  A  and  using  (B.26)  gives. 


dH{6) 


(N«/2) 


(h^(t  -  Ti)4>)‘ 


(m/2) 


h[t  -  ri)lT{t  -  Ti)  if) 


(iV„/2)^  ^ 

Taking  the  negative  of  the  expected  value  on  both  sides  gives. 


®\aA0A^J  (N„/2) 


(B.9) 
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Partial  differentiation  of  (B.5)  w.r.t  Ai  gives, 

“  Aih^{t  -  {-h^{t  -  n))  dt+ 


-  Ti)(t)){-Aih^ {t  -  Ti))  dt 


Taking  the  negative  of  the  expected  value  of  both  sides  and  using  (B.26)  gives, 

f  (9^/(0)  1  2Ai  r  ,Tu  ^  \J.uTu 


dAidff)^ 


\  =^/  h'^  {t  -  Ti)(f>h^  {t  -  Ti)  dt 

)  t\Q  Jq 


94  •  r'^ 

=  /  Kt  -  n)h^{t  -  Ti)  dt 

^^0  Jo 


(r/2)  ,1 

(N„/2) 


and  so 


^  f  Ti(e)  1  _  (r/2)  .  j 


(B^lO) 


Partial  differentiation  of  (B.5)  w.r.t  cf)  and  using  (B.26)  gives. 


dH{e) 

d4>d4>^ 


^'j'  M—l 


0  Jo 


{-Aih^{t  -  Ti))  {-Aih^{t  -  Ti))  dt 


{No/2) 


h(t  —  Ti)h7'{t  —  Ti)  dt 


and  so 


^  {T/2)  2 


^  f  dH{e)  \  (Til) ,  *49  ,,  (r/2)A^A, 

[d4>d4T]  (iV„/2)  <“-■>  h;  •  (]V„/2) 


(B.ll) 


Putting  (B.6),  (B.7),  (B.8),  (B.9),  (B.IO),  (B.ll)  back  in  (B.l),  we  have 


(2,rFo)2,^TLLT^(ji3^g(A))2  (2,rFo)(<#>2’L,^)(diag(A))  (2,rFo)(A  © 


(T /2') 

le  =  (2,rFo)(<#.^L^<^)(diag(A)) 

(27rFo)L'*’<#.(A  0  A)2’ 


(<^-‘  (^)Im 


(A  A)I(2Ar-i) 
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The  CRLB  matrix  for  the  unknown  parameter  vector  9  is  the  inverse  of  the  matrix 
Xq.  But  in  Appendix  B.4  it  is  shown  that  the  null  space  of  Xg  is  not  empty  and 
so  it  is  not  invertible.  This  is  because  the  log-likelihood  function  is  not  uniquely 
dehned  by  the  model  in  (23).  To  eliminate  the  over  parameterization  we  use  the 
following  transformations. 


T 

A' 

4>' 


=  [  (n  -  To)  (r2  -  To)  ■  •  •  (tm-I  - 

=  (lMo)[^l  Am-iY 


1 

0(l,2Af-2) 

0(2V-2,1) 

J-Af-l 

f-v-i 

Ivf-i 

— Ivf-i 

To) 


T 


diag(/i(-ro))(/) 


(B.13) 


Let  O'  =  [r'^  A'^  4>'^Y-  This  is  a  function  of  0.  Let  H 


dO' 


be  the 


Jacobian.  If  H  has  row  vectors  that  are  linear  combinations  of  those  eigenvectors 
of  Xq  that  have  nonzero  eigenvalues  then  the  CRLB  of  O'  is  given  by  [40]. 

The  t  is  used  to  represent  the  generalized  inverse.  This  condition  is  verihed  in 
Appendix  B.4.  Therefore, 


xY  =  hxJh^ 


(B.14) 


Alternately,  the  log-likelihood  function  for  this  model  with  the  transformed  pa¬ 
rameters  is  given  by 


rp  rp 

^^0  Jo  JVo  Jo  ^ 

(B,15) 

So,  computing  the  derivatives  to  hnd  the  FIM  as  done  previously  yields  the  FIM 
for  the  transformed  parameter  vector  as 


T/2 

(^Vo/2) 


(2,rFo)=^<#>' '^LL'^<#>'(diag(A'))^  (27rFo)(<#.' ^L,^')(diag(A'))  (27rFo) (A' 0  A')<^' ' 

(27rFo)(<#)' ^L^<#<')(diag(A'))  (<t>' <t>' )1m -1  A'<#>' 

(27rFo)L'r<#>'(A' O  A')'^  <l>' A' (1  +  A' A')I(2N-1)  . 


(B.16) 


We  have  verified  numerically  that  (B.14)  is  equivalent  to  (B.16).  The  elements  of 
the  TDOA  vector  t'  are  given  by 


,  I  .  a/ (x^  -  Xi)^  -h  iVj,  -  ViY  Yixj,  -xoY  -h  iVj,  -yoY 

A  —  (a  To)  — 

c  c 

so  that  the  new  parameter  vector  r'  is  a  function  of  only  the  emitter  location 
{x^,yY-  So,  if  we  let  rf'  =  [x^  y^Y  and  ex'  =  [r/'^  A'^  4>'Y^  have 


Xa’  — 


do' 

dex'  ^ 


X, 


e' 


dO' 

dex''^ 


dO' 

dex''^ 


dO' 
dex'  ^ 


(B.17) 
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The  Jacobian  is  given  by  the  (2M  +  2N  —  3,  M  +  2N)  matrix 


dO' 

da''^ 


dr' 

dr}' 

dM 

dr}' 

d(f)' 

dr}' 

dr' 

dr}' 


dr' 

dA’^ 

dA’ 

dA'^ 

dcf)' 

dA'^ 


dr' 

d<t)'^ 

dA! 

d4>''^ 

d(t>' 

d<p''^ 


LI(M,2V-1) 


0(Ar,2)  Im-1  0(m,2V-1) 

0(2V-1,2)  0(2Ar-l,M)  I(2V-1) 


1)  X  2  matrix  j 

= 

{x^  -  Xi) 

1 

0 

S 

1 

iS 

1 

1 

0 

di 

(x^  -  X2) 

{x^  Xo) 

di 

{Vt  -  J2) 

{Vt  -  Vo) 

d2 

do 

(J2 

do 

Xrp  1) 

(Xj,  -  Xo) 

{x^  —  Xm-i) 

'5' 

1 

0 

(B.18) 


(B.19) 


L  dM-i  do  dM-i  do  J 

where  di,  i  =  0,  •  •  •  ,  M  —  1  is  the  distance  between  the  sensor  i  and  the  emitter. 


B.1.2  Signal  known  with  unknown  transmission  time 

From  (34)  we  have  the  log-likelihood  function  as 


1  pT  M-l 

~W  /  “  Aih^{t  -  dt 

^''0  Jo 


KC)  =  -^ 


(B.20) 


Here  the  2M  x  1  unknown  parameter  vector  is  ^  =  [r^  A^]^.  So,  the  FIM  is  given 
by 

r  p  r  \  i 

Xdrdr^]  {drdA^j 

=  (B.21) 

p  r  d\o)  \  p  r  dH{e)  \ 

^  1  dAdr^  J  ^  1  dAdA^  J 
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Using  (B.6),  (B.7)  and  (B.9),  we  have 


(T/2)  '  (27ri^o)^0  LL  0(diag(A))2  (27rFo)(0^L(/))(diag(A))  ' 

~  i  M  (B.22) 

[  (27rFo)(0^L'^0)(diag(A))  (0'^0)Im 


B.2  Maximum  Likelihood  Estimator 

Here  we  will  derive  the  MLE  for  the  two  cases  of  signal  unknown  with  unknown 
transmission  time  and  signal  known  with  unknown  transmission  time. 

B.2.1  Signal  unknown  with  unknown  transmission  time 

The  log-likelihood  function  with  the  transformed  parameters  is  given  by 

rjn  rp 

^^0  Jo  JVo  Jq  ^ 


Partial  differentiation  w.r.t  cf)'  gives, 


(B.23) 


^ojo  — 


2  (u(t)  -  A[h^{t  -  {-A[h^{t  -  t'))  dt 


In  order  to  hnd  the  maximum,  we  equate  the  above  partial  derivative  to  zero, 
which  gives. 


ro{t)h^  (t)  dt 


ri{t)h^ {t  —  t[)  dt 


(/  ~  d^  =  0 


Using  the  properties  of  the  vector  h{t)  as  shown  in  Appendix  B.3,  we  have 
!o  ^o{t)h^{t)  dt  -  (T/2)0'^+ 

E"r"  -  7i)  dt)  -  (T/2)  E"7'  -4' =  0 

If  we  replace  the  integrals  with 

Yq  =  f  ro{t)h{t)  dt  and  y'i  =  f  ri{t)h[t  —  t[)  dt  for  i  =  1,  2,  •  •  •  M  —  1, 

Jo  Jo 

then  we  have 


y'/  -  {T/2)d>''^  +  5^  4'y'^  -  (T/2)  =  0 
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So,  the  MLE  of  (p'  is 


(2/T)  (y'o  +  E£7' 4y0 

4>  =  — 


(l  +  A'^A'^ 


Putting  this  back  in  (B.23),  we  have 


M-l 


xl{t)  dt^- 


1  (2/T)  (yi^  +  E,"r'  (y'„  +  ESt'  -4'y' 


A^n 


(1  + A'^A'" 


Maximizing  /(0^)  w.r.t  A'  and  r'  is  equivalent  to  maximizing  the  second  term.  So, 
let 


/(A',t')  = 


M-l 


+  E:=i 


A'yE)(yo  +  E£7^A'y' 


(l  +  A'^A' 


If  we  let  Y'  =  [yg  y'^  ■  ■  ■  yV-i]  be  the  2 A  —  1  x  M  matrix  then  the  maximum  value 
of  /(A',  r')  w.r.t  A'  is  fma.xi'’’')  is  equal  to  the  maximum  eigenvalue  of  Y'Y'  Let 
B'  =  YY'^.  The  matrix  B'  is  a  function  of  r'  =  which  is  a  function  of 

the  emitter  location  {x^,y^).  So,  the  MLE  of  the  emitter  location  is  found  by 
maximizing  the  maximum  eigenvalue  of  i.e. 


{xT^yx)  =  argmaxAmax(B'(a;j,,|/j,))  (B.24) 

(xT,yT) 


B.2.2  Signal  known  with  unknown  transmission  time 

The  log-likelihood  function  is  given  by 

1  nX  M-l 

/(C)  =  “  Aih^{t  -  dt  (B.20) 

AoJo 

Partial  differentiation  w.r.t  Ak  gives, 

[  2  (rfc(t)  -  Akh^{t  -  Tk)4>)  -  Tk)4>)  dt  =  0 

for  each  of  A;  =  0, 1,  •  •  •  ,  M  —  1.  Equating  this  to  zero  to  hnd  the  maximum  value 
gives, 

(y  rk{t)h{t  -  Tk)  d^  -  AkcfF  h{t  -  Tk)h^ {t  -  Tk)  di^  (p  =  Q 

If  we  replace  the  integral  with  y^  =  rk{t)h(t  —  Xk)  dt  and  use  the  properties  of 
the  vector  h{t)  as  shown  in  Appendix  B.3,  we  have 

cp^jk  -  {T/2)Ak(P^(P  =  0 
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So,  the  MLE  of  is 


Ak 


(T/2)0^(/>’ 


k 


-  ,M  -  1 


Putting  this  back  in  equation  (B.20),  we  have 


^(0 


N, 


M—l 

/  rf(t)  dt 
»  i=o  •'» 


2 f  •ry. 

\(TI2)4?'4, 

2 


ri{t)h^ it  -  Ti)  dt]  <p+ 


) 


N, 


M-l 

E 


M-l  /  ,T  T  ,  \ 


0  .^0  Jo 


ivi-L  /  rr  T  / 

rf(«)  d*  +  -  y  f 


Maximizing  /(^)  w.r.t  r  is  equivalent  to  maximizing  the  second  term.  So  the  MLE 
for  1}  is  given  by 

M-l 

f)  =  argmaxN  =  argmax0^B0  (B.25) 

^  t"o  ^ 

where  B  =  YY"^  and  Y  =  [yg  y^  ■  ■  •  yM-i]  with  y^  =  ri{t)h{t  -  n)  dt,  i  = 
0, 1,  •  •  •  M  —  1.  B  is  a  function  of  {x^,y^,to). 


B.3  Properties  of  h{t) 

The  time  dependent  vector  h{t)  that  was  used  for  modeling  the  problem  in 
equation  (21)  has  some  interesting  properties  which  simplihes  the  derivation  of  the 
CRLB  and  the  MLE.  These  properties  are  derived  here.  We  have 


hit  -  Ti) 


cos27rFo(^  ~  t)  cos27r(Y  —  l)Fo(^  —  p)  sin  27rFo(^  ~  p) 
■  ■  ■  sin  2t{N  -  l)Fo(^  -  n)f 


Differentiating  both  sides  w.r.t  Ti,  we  get 


dh(t  —  Ti) 


dri 


=  2t:Fq  [  0  sin  27rFo(t  —  Ti)  2  sin  27r2Fo  (t  —  Tj)  •  •  •  (iV  —  1)  sin  2't{N  —  l)Fo(t  —  Ti) 

—  cos  27rFo(t  —  Ti)  —2  cos  27r2Fo  (t  —  r^)  •  •  •  —  (W  —  1)  cos  27r(iV  —  l)Fo(£  —  r^)  ] 


Let 


0 


{N,N) 


0(i,Ar-i) 

diag(l,2,  •••  ,  Y  -  1) 


—  [  0(Ar-i,i)  diag(l,  2,  •  •  •  ,  Y  —  1)  ]  0(Ar_i^Ar_i) 

So,  we  have  the  partial  derivative  of  h{t  —  Ti)  w.r.t  to  Ti  as 

dh{t  -  Ti) 


dTi 


=  (27rFo)Lh,(t  -  Ti) 
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Next  we  will  compute  the  integral  h{t)h^ (t)  dt.  We  have, 


h{t)h^  {t)  = 


r  1  1 

V2 

r  1  1 

V2 

COS  27rFot 

COS  27rFot 

cos27r(N  —  l)Fot 

cos  2Tr{N  —  l)Fot 

sin  27rFot 

sin  27rFot 

sin27r(N  —  l)Fot 

sin27r(N  —  l)Fot 

T 


Let  us  compute  each  of  the  integrals  in  this  2N  —  lx  2N  —  1  product  matrix 
separately.  Integral  of  the  hrst  element  is 

Integrals  of  the  elements  on  the  diagonal  are  given  by 


and 


cos^  2nkFot  dt  =  - 


(1  +  cosAirkFot)  dt 

\T 


t  + 


sin  AirkFot 
ATlkFn 


sin^  2TTkFQt  dt  =  - 


2 

T 


(1  —  cosAirkFot)  dt 
sin  AnkF^F 


AirkFo 


0 


for  /c  =  1,  2,  •  •  •  iV  —  1.  Integrals  of  the  rest  of  the  elements  are  given  by 


/„  cos  27r/cFof  sin  27rnFot  dt  =  -  sin  27r(n  +  k)FQt  dt  +  sin  27r(n  —  k)FQt  dt 


1 
2 
_  1 
“  2 
=  0 


cos  27r(n  +  k)FQt  cos  27r(n  —  k)FQt 
47r(n  +  k)FQ  47r(n  —  k)FQ 


■^T 


and 

rT 


Jq  cos  2nkFQt  cos  2nnFQt  dt  =  - 


1 
2 
_  1 
“  2 
=  0 


cos  27r(n  +  k)FQt  dt  +  cos  27r(n  —  k)FQt  dt 

T 


sin27r(n  +  k)FQt  ^  sin27r(n  —  k)FQt 


47r(n  +  k)FQ 


47r(n  —  /c)Fo 


0 
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and 


T  / 

Jp  sin  27ikFQt  sin  27rnFot  dt  =7:  cos  27r{k  —  n)Fot  dt  —  cos  27r(n  +  k)FQt  dt 

^  -Jo 

1  sin27r(/c  —  n)Fot  sin27r(n  +  A;)Fot 

2  47r(/c  —  n)Fo  47r(n  +  k)FQ  ^ 

=  0 


for  A;,  n  =  1,  2,  •  •  •  ,  iV  —  1.  Therefore,  we  have  the  integral  of  h{t)h7' {t)  as  a  scaled 
identity  matrix  given  by 


1  0 

rT  T  0  1 

h{t)h^{t)  dt  =  — 


-  {T/2)I(2N-i) 


[0  0  ■■■  IJ 

Since  h(t)  is  periodic  with  period  T,  for  any  r, 

[  h{t  -  T)h^{t  -T)dt=  [  h{t)h^{t)  dt  =  (T/2)I(2jv-i)  (B.26) 

Jo  Jo 

Now,  we  compute  the  integral  of  the  cross-product  of  the  partial  derivatives  of 
h(t  —  Ti)  w.r.t  to  Ti 

dh{t  -  Ti)  dh^  it  -  Ti) 


=  /q  (27rFoLAi(t  -  r^))  {27rFoLh{t  -  Ti))  dt 


(27rFo)^L  /J’  h{t  -  Ti)h^ {t  -  n) 


i)  dt  \  L 


=  (T/2)(27rFo)2LL^ 

(B.27) 

and  the  integral  of  the  cross-product  of  h(t  —  Ti)  with  its  partial  derivative  w.r.t 
to  Ti  is 

[  — —h^{t  -  Ti)  =  27rFoLh{t  -  Ti)h^ {t  -  Ti)  dt 


(27rFo)L  h{t  -  Ti)h^ {t  -  Ti) 


(B.28) 


=  (T/2)(27rFo)L 


B.4  Transformation  of  the  Parameters 

In  Section  3.3  we  discussed  the  relationship  between  the  unknown  attenuation 
factors  and  the  unknown  signal,  and  between  the  unknown  TOAs  and  the  unknown 
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signal.  Here  we  will  show  that  the  FIM  given  in  equation  (24)  is  rank  two  dehcient. 
Then  we  will  show  that  the  transformation  given  in  equation  (25)  satishes  the 
conditions  given  in  [40].  From  (24),  we  have 

(2,rFo)2<^'^LL'^<#>(diag(A))2  (2,rFo)(<^^L,^)(diag(A))  (27rFo ) (A  ©  A) ,^'^L  ' 

(27rFo)(<#)^L^<^)(diag(A))  A0^ 

(27rFo)L'^</.(A0  A)^  ^A'^  (A^A)I(2„_i) 

If  i/i  =  [  llj  -(0+(27rFo)L'^0)  f  and  1/2  =  [  -A^  (0- (27rFo)L^(/)) 

then  it  can  be  verihed  that  XgUi  =  0  and  Xgi/2  =  0.  Therefore  Ui  and  U2  are  in  the 
null  space  of  Xg.  Also,  Xg  +  (l/2)i/ii/f  +  {1/2)1/21^2  is  non-singular.  This  means 
that  ui  and  1/2  are  the  basis  vectors  for  the  null  space  of  Xg  and  so  the  matrix  Xg 
is  rank  two  dehcient.  The  Jacobian  of  the  transformation  is  given  by. 


(T/2) 

(A'o/2) 


H  = 


-  dr' 

dr' 

dr'  - 

'M 

(90 

dA' 

dA' 

dA' 

dr 

dA 

(90 

d(t)' 

d(t)' 

(90' 

L  dr 

dA 

(90  . 

(B.29) 


Now,  we  will  compute  each  of  the  derivatives  in  the  Jacobian  matrix.  The  elements 
of  the  hrst  sub-column  are  given  by. 


dr' 

dr' 

(90 


—  ^  ([— Ivf-i  Im-i]t) 
=  [— 1m-i  Im-i] 

=  0i^m-i,m) 

=  0(m_i,27V-1) 


The  elements  of  the  second  sub-column  are  given  by. 


(9A' 

dA' 


dA' 

(90 


=  ^  (([0(M-i,i)  Im-i]A)  (ef a) 

=  (i-Mo)[~A'  Im-i] 

=  0(M-1,2V-1) 
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where  ei  is  the  first  column  of  an  M  x  M  identity  matrix.  The  elements  of  the 
third  sub-column  are  given  by, 


d<p' 

di- 


d<p' 

dA 


d<p' 

d<p 


d 

1 

0(l,2Ar-2) 

—  Ao 

0(2Ar-2,l) 

1  1 

I77-I 
—  I77-I 

][ 

1 

0(1,277-2) 

- 

Aq 

0(277-2,1) 

1  1 

Ia7-1 
—  Ia7-1 

1 

0(1,277-2) 

- 

Aq 

0(277-2,1) 

1  1 

Ia7-1 
—  Ia7-1 

diag(h.(— tq))^ 


dh{  —  TQ)  \ 

diag  (  - - j  <j)  0{2Ar-l,l)  ■■■  0{2Ar-i,i) 


(27rFo)Ldiag(h.(-To))<^  0(2N'-1,1)  •••  0(2N'-1,1) 


dA 


(e^A) 


1  0(l,2iV-2) 

n  r  IiV-1  IaT-I 

0(2N-2,1)  [ 


diag(h,(-To))0 


0 


dtp 


^(l,2N-2) 

,  liv-1  Ia^-i 
(2Ar-2,l)  [ 

1 


(e^A) 


0(l,2Ar-2) 

n  r  ^N-1  ^N-1 

0(2iV-2.1) 


diag(h,(-To))0  0(2Ar_i,i) 
diag(fi(-To))0^ 


-J(2Ar-l,l) 


^0 


^(2Ar-2,l) 


^{l,2N-2) 
IaT-I  lAf-1 
Iat-i  — Iat-i 


diag(Ai(-To))  0(2Ar-l,l) 


^{2N-1,1) 


For  the  transformed  parameters  to  have  hnite  variance,  the  row  vectors  of  H  must 
be  equal  to  the  linear  combinations  of  those  eigenvectors  of  Xg  that  have  nonzero 
eigenvalues  [40] .  In  order  to  show  that  the  row  vectors  of  H  are  linear  combinations 
of  those  eigenvectors  of  Xg  that  have  nonzero  eigenvalues,  it  is  enough  to  show  that 
the  row  vectors  of  H  are  orthogonal  to  the  null  space  of  Xg.  That  is,  it  is  enough 
to  show  that  Hi^i  =  0  and  Hi/2  =  0.  Now, 


Hi/i  = 


dA' t  dA'  A 


^{<1,  +  (2rf„)L^</.) 


dp' 

dT 


1m  + 


Substituting  the  partial  derivatives  that  we  computed  previously  and  further  sim¬ 
plifying  gives. 


[— liU-l  I(Ar-l)]lM  +  0(M_i4)  —  0(M-1,1) 

0(M-i,i)  +  (lMo)[— A'  I(a^_i)]A  —  0(A^_i^i) 

24oP(27rFo)Ldiag(/i(-ro))0  -h  24oPdiag(/i(-ro))0- 

(24oPdiag(h,(-ro))0  -h  24oPdiag(h,(-ro))(27rFo)L'^(/)) 


where 


1 

0(l,2Af-2) 

P  = 

0(27V-2,1) 

I(W-l)  I(Af-l) 

I(W-l)  — I(V-l) 

Using  the  fact  that  Ldiag(h,(— tq))  =  diag(/i(— ro))L^,  we  have  Hi/i  =  0.  Similarly 
it  can  be  shown  that  iii'2  =  0. 


74 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED. 


C  Optimal  Sensor  Configuration 


Here  we  will  analyze  the  effect  of  the  sensor-emitter  geometry  on  the  FIM  for 
the  emitter  location.  We  will  derive  the  D-optimal  configurations  for  the  uncon¬ 
strained  and  the  constrained  cases  that  are  described  in  Section  4.1.  From  (30) 
and  (B.18)  we  have  the  FIM  for  the  emitter  location  vector  r]'  =  [a;^  y^]  as 


X^/  — 


dr' 
dr)'  ^ 


(27rXo)V'^LL^0'(diag(A'))^ 


(C.l) 


If  we  assume  that  the  signal  level  is  the  same  at  all  sensors,  i.e,  Ai  =  A  for  all 
then  we  have  diag(A')  =  AI(m-i),  where  I(m-i)  is  the  (M  —  1)  x  (M  —  1)  identity 
matrix,  and  we  have  the  above  FIM  as 

I,,  =  (2rf„)V' (C,2) 


If  we  let  tfji  be  the  angle  of  sensor  i  from  the  positive  x-axis  measured  at  the  emitter 
as  shown  in  Figure  (D.l),  then  we  can  write  (B.19)  as 


)  =  (l/c) 


drf 


(cos'll  —  COS'lpo) 
(cos '^2  —  COS'^o) 

{cOS-lpM-l  -  COS'lpo) 


(sin'll  —  sin'^o) 
{simjj2  -  simjjo) 

{sin'ipM-i  -  sm-ipo) 


Now  let 


(C.3) 


Figure  C.l.  Definition  of  the  angle  ipi. 


J 


cos  '^0 

sin'^0 

1 


cos'll 

sin'll 

1 


cos  iJM-i 
sin'ipM-i 
1 
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Using  elementary  operations,  we  can  transform  J  to 


COS'^0 

cos  ^jJl  —  cos  '^0  •  • 

■  COS'lpM-l- 

cos  fjQ 

f  dr'  V  ■ 

siii'^o 

sin  '^1  —  sin  ipQ  ■  ■ 

■  smfjM-i- 

sin  fjQ 

= 

U  C  1 

\dr]''^  ) 

1 

0 

0 

_  1 

0^ 

where  u  =  [cos'^o  sin'^o  1]^  and  0  is  the  (M  —  1)  x  1  zero  vector.  Now,  since 
the  elementary  operations  do  not  change  the  determinant,  we  have 


det(Jj'^)  =  det 


u  c 


dr' 

dr^'T 


T  -\ 


U" 

Ot' 
dr]'  ^ 


1 

0 


=  det 


uu^  + 


dr' 

dr]''^ 

T 

U 


dr' 
dr]'  ^ 


u 

1 


=  det 


+  c^ 


So,  det 


(1/c^)^  det(Jj'^)  so  that  we  have. 


det(X,,)  =  hWWdet(Jjy 


Maximizing  the  determinant  of  the  FIM  is  equivalent  to  maximizing  det(JJ^ 


C.l  Unconstrained  Geometry 

This  problem  can  be  stated  as  follows: 
We  have  g  :  — )■  such  that 


g{xl))  =  det(JJ^) 

where  V’  =  [t/’o  t/’i  •  •  •  'h  =  [— tt,  tt]^  is  a  convex  subset  of  .  Find 


arg  max  g{'ip) 


We  have. 


JJ^  = 


EM—1  2  /  I  •  /  I 

cos^Ui  2^i=o  cosUiSmUj  cosf, 


^M-l 


M-1 


EM—1  I  •  I  •  2  /  •  / 

.^0  cos'diSim/ji  sm  ^/Ji  sim/j, 


M-l  ■  2 


.M-1  . 


V— .M— 1  / 

Ei=o  COSl/>i 


Ei=o  sin^i 


M 
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By  Hadamard’s  inequality, 


(M-l  \  /M-l  \ 

cos^  ipij  f  sin^  ipij  (M) 

The  equality  holds  if  and  only  if  JJ^  is  diagonal.  That  is  iff 

cos  ipi  sin  ipi  =  0 

E£ycos^,  =  0  (C.4) 

E*=o^sin^i  =  0 

Also, 

(E*=o  ^  cos^  V'i)  (Ei^o  ^  sin^  =  {YZYih  +  2Cos2V’i))  (E*=o^(5  -  5  cos  2^’*)) 

=  {f  + I  YZY  cos2V’i)  (f  -  5  YZY  cos2V’i) 

=  ¥-i(E.=Ecos2V-.)'<^ 

This  inequality  holds  iff 

M-l 

y^cos2-^j  =  0  (C.5) 

i=0 

When  the  sensors  are  placed  around  the  emitter  at  equi-angular  distances  (as 
shown  in  Figure  17(a)),  then  we  have  =  27r  ,i  =  0,  •  •  •  ,M  —  1.  This 

conhguration  satishes  the  conditions  in  (C.4)  and  (C.5).  Therefore,  this  is  an 
optimal  conhguration.  Also,  when  M  is  a  multiple  of  three,  if  the  sensors  are 
distributed  equally  at  the  three  vertices  of  an  equilateral  triangle  inscribed  in  the 
circle  around  the  emitter,  as  shown  in  Figure  17(b),  then  this  conhguration  also 
satisfy  these  conditions  and  hence,  is  also  an  optimal  conhguration. 

C.2  Constrained  Geometry 

This  problem  can  be  stated  as  follows: 

ITe  have  g  :  — )■  such  that 

gi'il))  =  det(JJ^) 

where  -0  =  [t/’o  t/’i  •  •  •  f/’M-i]^?  d'  =  [—6,6]^  is  a  convex  subset  of  .  Find 

arg  max  gl'il)) 

Notice  that  here  the  domain  is  constrained.  For  this  constrained  geometry  problem 
we  will  derive  the  D-optimal  conhgurations  for  three  and  four  sensors. 
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Three  sensor  case 

Here  we  use  a  simple  transformation  of  variables.  Let  cosipi  =  Xi,  i  =  0, 1,  2, 
so  that  sin'll  =  a/1  —  xf.  With  this  transformation  of  variables,  the  problem  can 
be  restated  as: 

Suppose  that  x  =  [xq  Xi  Xq\  and  X  =  C]  x  [xq,  C]  x  C]  C  Let 


J  = 


Xq  Xi  X2 


f{Xo)  f{xi)  f{x2) 
111 


where  f{x)  =  y/1  —  x‘^  and  let  g  :  ^  ^  such  that  g{x)  =  det(JJ^)  =  (det(J))^. 

Find 

arg  max5f(a;) 
xex 


The  setup  is  as  shown  in  Figure  C.2. 


Figure  C.2.  Three  sensor  setup. 

Now,  for  arbitrary  Xi  >  we  have  Xq  G  [— C,a;i].  We  will  show  that  for  all 
arbitrary  hxed  Xi,  X2  G  X,  g{xo)  is  a  convex  fnnction  of  xq  for  xo  G  [— C,a;i].  We 
have 

Now, 


det(J)  =  xoifixi)  -  f{x2))  -  f{xo){xi  -  X2)  +  {,xif{x2)  -  X2f{xi)) 


d 

dxQ 


det(J)  =  {f{xi)  -  f{x2))  -  fixo){xi  -  X2) 


dgjxo) 

dxo 


2det(J)((/(a;i)  -  f{x2))  -  f{xo){xi  -  X2)) 
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^  =  2  +  2det(J)?!^ 


Now 


Oxq 


(9det(J) 


dxQ 


dxn 


=  ifixi)  -  fix2))  -  fixo)ixi  -  X2) 


<92  det(J) 


dxl 


=  -f"{xo){xi  -X2)  <0 


since  /  is  concave  and  xi  <  X2.  Next  we  will  show  that  det(J)  <  0 

det(J)  =  xo{f{xi)  -  f{x2))  -  f{xo){xi  -  X2)  +  {xif{x2)  -  X2f{xi)) 

=  f{xo){x2  -  Xi)  +  f{x2){xi  -  Xq)  -  f{xi){x2  -  Xq) 

Choose  a,  (3  &  (0,1)  such  that  a  +  (3  =  1  and  axo  +  /3x2  =  xi.  So,  we  get 
a  =  (x2  —  Xi)/{x2  —  xo)  and  /3  =  (xi  —  xo)/(x2  —  xq).  Now,  since  /  is  concave, 

/(..)  =  /(a.„  +  /?x.)  >  «/(x„)  +  W(x.)  = 


(X2  -  Xo) 


(X2  -  Xo) 


(X2  -  Xo)f(Xi)  >  f(xo)(x2  -  Xi)  +  f(x2)(xi  -  Xo) 
d^ff(xo) 

Therefore,  det(J)  <  0  and  so  — —  >  0.  This  means  that  g(xo)  is  a  convex 

OXq 

function  of  xo  G  [— C,a;i].  The  maximum  of  g{xo)  occurs  at  xo  =  —(  or  xo  =  xi. 
When  Xo  =  xi,  g{x.)  =  0.  So  the  xo  that  maximizes  gigi)  is  xo  =  — C-  Similarly,  we 
can  show  that  the  X2  that  maximizes  ^'(x)  is  X2  =  (■  Now  £x  xo  =  — C  and  X2  =  (■ 

So  xi  e  [-C,C]- 

Now, 


det(J)  =  Xoifixi)  -  f{x2))  -  f{xo){xi  -  X2)  +  {xif{x2)  -  X2f{xi)) 


det(J)  =  ifix2)  -  fixo))  +  f'ixi){xo  -  X2) 

UXi 

=  2det(J)((/(a;2)  -  f{xo))  +  f{xi){xo  -  X2)) 

UXi 

Equating  this  to  zero  to  determine  the  extrema 

2det(J)((/(a;2)  -  f{xo))  +  f\xi){xo  -  X2))  =  0 


det(J)  =  0  or  {f{x2)-f{xo))  +  f{xi){xo-X2)  =  0 
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If  det(J)  =  0  then  g^-x)  =  0,  so  this  is  not  a  maximum.  When  {f{x2)  —  /(a^o))  + 
f'{xi){xo  —  X2)  =  0,  we  have 


f'M 


{f{x2)  -  f{Xo)) 

{X2  -  Xo) 


/(C)  -  /(-C) 
C  +  C 

Therefore  xi  =  0.  The  solution  for  the  problem  is 


argmax5((x)  =  {[-(  0  C]'^} 
xex 


Four  sensor  case 

Here  we  assume  the  domain  T  =  [0,  tt]  x  [ipQ,  tt]  x  [/’i,  tt]  x  [ip2,  tt]  C  3?"^.  Notice 
that  here  we  are  assuming  that  the  constrained  arc  is  a  semicircle.  The  more 
general  case  where  0  <  ipo  <  ipi  <  %Ij2  <  n  —  6  where  0  <  6^  <  7r/2  appears  to 
be  more  difficult.  Using  the  argument  from  above,  we  can  show  that  for  optimal 
conhguration  =  0  and  '^3  =  tt.  This  is  shown  in  the  Figure  C.3.  So,  we  have 


J  = 


0  sin'll  sin '^2  0 

11  11 


We  want  to  prove  that  (1/4)  det(JJ^)  <  2.  Note  that  from  Cauchy-Binet  theorem 


det(Jj'^)  =det' 


det^ 


I  cos  t/ji  cos  '<p2 
0  sin  sin  '^2 

II  1 
1  cos '^2  —1 

0  sin  'ip2  0 
1  1  1 


det^ 


1  cos  ipi  ■ 

det^  (  I  0  sin'll  0 
111 
cos'll  cos '^2  —  1 

sin'll  sin '^2  0 

1  1  1 
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=  (sin  ipi  +  sin  '^2  cos  ipi  —  sin  tjji  cos  ifj2  —  sin  ■^2)^  +  (2  sin  ipiY  +  (2  sin  '^2)^+ 

(sin  'ip2  cos  ^jJl  —  sin  ■^i  +  sin  ip2  —  sin  ipi  cos  '^2)^ 
=  4sin^  ipi  +  4sin^  '^2  +  2(sin'^i  —  sin ■^2)^  +  2(sin'^2  cos'll  —  sin'll  cos '^2)^ 

Hence  we  need  to  prove  that 

sin^  V’l  +sin^  V’2  +  (l/2)(sin  V’l  —  sin  V'2)^  +  (l/2)(sin  ^^2  cos  V’l  —  sin ■0i  cos  V'2)^  <  2  (C.6) 


(C.6)  is  equivalent  to  the  following  inequality,  which  can  be  easily  verihed 


(1/2)  (sin  ^/)2  cos  V’2+sin^/i  cos  V’i)^  +  (l/2)(sin^  V’l+sin^  ■02  —  2)^  >  (1/2)  (sin  0i 
Now,  if  we  show  that 


sin  02)^ 
(C.7) 


(sin^0i  +  sin^02  —  2)^  >  (sin0i  —  sin02)^  (C.8) 

then  (C.7)  is  proved  since  (sin02cos02  +  sin0icos0i)^  >  0.  Without  loss  of 
generality  (  due  to  symmetry  of  (C.8)  ),  we  assume  that  sin  02  >  sin0i.  Hence 
(C.8)  is  equivalent  to 


2  —  sin^  01  —  sin^  02  >  sin 02  —  sin  0i  (C.9) 

which  is  equivalent  to 

5/2  -  (sin 01  -  1/2)2  _  ^  ^^2)^  >  0  (C.IO) 

Since  (sin0i  —  1/2)^  <  1/4  and  (sin 02  +  1/2)^  <  9/4,  inequality  (C.IO)  is  easily 
verihed  and  this  proves  inequality  (C.6).  Equality  in  (C.IO)  holds  when  sin0i  =  0 
or  1,  sin 02  =  1  under  which  condition  the  equality  in  (C.7)  holds.  This  means  that 
sin  01  =  0,  sin 02  =  1  and  sin0i  =  1,  sin  02  =  1  are  optimal  solutions.  Note  that 
since  we  have  assumed  sin  02  >  sin0i,  we  would  also  have  sin0i  =  1,  sin  02  =  0 
as  the  optimal  solution.  Therefore  the  optimal  sensor  conhgurations  are  given  by 

argmax5((0)  =  { [0  0  7r/2  tt]^,  [0  7r/2  7r/2  tt]^,  [0  7r/2  tt 
061' 
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D  Derivation  of  the  Fisher  information  matrix 


Assume  a  stationary  emitter  and  stationary  receivers.  Also  to  simplify  the 
derivation,  assume  a  known  signal  and  known  transmission  time.  We  wish  to  de¬ 
termine  the  CRLB  for  the  emitter  position. 


At  receiver  i,  [i  =  1,  -  ,  M)  we  receive 

ri(t)  =  Ai{x^,y^)s{t  -  Ti{x^,y^))  +  Wi(t) 

where  s(t)  is  known,  Aj (x j, ,  )  and  Ti{Xj,,yj,)  are  known  functions  of  {x^,y^)  and 

Wi(t)  is  white  Gaussian  noise  with  the  same  power  spectral  density 

Pw{F)  =  Nq/2  —  oo  <  F  <  oo 


at  each  receiver.  The  noises  are  independent  from  receiver  to  receiver.  We  observe 
r(t)  =  [ri(t)  ■  ■  ■  rM(t)],  0  <t  <T.  Then  the  probability  density  function  is 

p{r{t);x^,y^)  =  coexp  ((-l/iVo)(5) 

where  Cq  is  a  constant  and 


M 

Q  =  '^  (ri(t)  -  Ais{t  -  Ti)Y  dt 
i=i  Jo 

and  where  Aj  and  r*  both  depend  on  the  emitter  position.  The  log-likelihood 
function  is 

L{x^,y^)  =  {-l/No)Q{x^,y^) 

To  hnd  the  CRLB,  the  FIM  matrix  is 


Now, 

dQ 

dx^ 


X{x^,y^)  =  {1/Nq)E 


d^Q 

d^Q 

dx"^ 

T 

dx^dy^ 

d^Q 

d^Q 

dx^dy^ 

dyl 

Ais{t 


t)) 


dsjt  -  Tj) 
dx^ 


dt 


d^Q 

dx‘^ 

T 


=  -2 


M 

E 

i=l 


f'T  Q 

{ri(t)-Ais{t-Ti))  — 

^  UXri 


-A, 


ds{t  -  Ti)  dAi 


dxr. 


dxr, 


S(t  -  Ti) 


dt-\- 


82 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED. 


i=i 


After  taking  the  expected  value,  the  hrst  integral  is  zero.  Thus, 


^  / d‘^Q\  ^  r  .  ds{t  —  Ti)  dAi 

r  1 

ox^  ox^ 

The  other  terms  are  found  similarly  so  that 

^  r  lo  9u  i^T  ^Vt)  dt  lo  aS  i^T  ,yT)  dt' 

Iix,,y,)  =  i2/No)J2 

L  Jo  aS {Xt ^Vt)  dt  /(f  9^22  {Xt ^Vt)  dt  . 

where  for  example 

ii)r  ^  ds{t-Ti)  dAi  ]  [  ds{t-Ti)  dAi 

a2l{XT:yT)  =  - h7^s(t-ri)  Ai — - h  7^ — s{t-Ti^ 


Converting  to  the  frequency  domain  by  using 


g{t)h{t)  dt  =  /  G*{F)H{F)  dF 


We  have  upon  letting  a;^  =  rji  and  =  rj2  so  that 


/o  g)Jn  dt  =  fJl  (^AiF  I  I  +  ^F  {s(t  -n)}) 

But  .7{s(«-n)}  =  S(F)exp(-j2rfTj)  and.7|  |  =  S(F)jf-exp(-j2irFTi) 


tTjmn*  =/!t,  tliS*(F)5|;;exp(j2xFTi)  +  ^S"(F)exp(j2xFri)  . 

A.S{F)j^  exp(-i2xFTi)  +  |^S(F)  exp(-j2xFT.)]  dF 

=  |S(f’)P  (2li|g;(j2xF)exp(j2xFT.)  +  |!jexp(j2xFTj))  . 

{Ai§^{-j27TF)exi>(-j2nFTt)  +  ||j  exp(-j2xFTi))  dF 

=  n  |S(f’)P  (A(j2xF)g;  +  1^)  (a.(-j2xF)|!1  +  |A)  dF 


=  IZ \S{F)?A!(2^Ff^lg  dF  +  |S(F)p|i|A  dF 


drjm  drjn 
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The  cross  terms  are  zero  since 


F\S{F)\^  dF  =  0 


So, 


where 


f  ,  f°°  .  „x2ir</r^\l2  it^a2^F  dTi  ^dAidAi 

/  gl2^dt=  /  {2T[Ff‘\S{F)\^  dF  Aj- — - V  S- — - — 

/o  7-00  dgmdgn 


S=  /  \S{F)\^  dF 


is  the  energy  of  the  signal. 

rT 


where 


ro  ,  a9  dTi  dTi  ^  dAi  dAi 

dt  =  SF^A^,^  ^  +S^  \ 

dr]m  drjn  drj^  drjn 

jrj27rFf\SiF)\UF 

r  = 


jTjS(F)\^dF 

is  the  mean-square  bandwidth.  But 

Ti  =  {l/c)^J{x^  -  XiY  +  {y^  -  ViY 


dTj 

dXrr, 


=  (1/c) 


{x^  -  Xi) 


a/ (Xj,  -  XiY  -h  iVj,  -  Vif 


=  (1/c)  COS^/>i 


dTi 


and  TT —  =  (1/c)  sin'll  where  ijji  is  the  angle  of  sensor  i  from  the  positive  x-axis 
OVt 

measured  at  the  emitter  as  shown  in  Figure  D.l.  Hence 


Figure  D.l.  Dehnition  of  the  angle  'di- 


[X(Xj,,  Hrp)  j 


1 


E".  -1  FF^AUl/c^ 


(Ko/2) 


COS^  Tpi 

sin  'ipi  cos  tfji 


sin  Tpi  cos  Tpi 
sin^  'ipi 


+ 


C  dAj  dAj 
^  drim  drin 
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or  finally  we  have 


X^Xt^Vt) 


£ 

(Wo/2)  c2 


E;U^n=:r,Vr) 


cos^  'ipi  sin  -ipi  cos  'ipi 
sin  'j/Zj  cos  ipi  sin^  ■j/'i 


+ 


£  s^M 
(Wo/2)  Z^i=l 


(D.l) 


D.l  With  no  Obstacles 

Assume  direct  path  to  all  sensors  and  no  multipath.  So,  the  A^’s,  the  attenu¬ 
ation  factors  can  be  modeled  as  inversely  proportional  to  the  distance,  i.e, 


-  xA  +  (vt  -  ViY  A 


where  Gi  is  a  constant. 


dx^ 

and  similarly, 


dXrr 


A  Vix^-xA  +  iy^ 


dAi{x^,yA  _  -Gj 


Putting  these  in  (D.l),  we  have 
AXt^Vt)  =(wf7^E*=l 


dyj 


sin  ipi 


f  f2  Gf 

COS^  -Ipi 

sin  ijji  cos  ipi 

|c2-R2 

sin  -ipi  cos  'ipi 

sin^  ipi 

+ 


COS  tpi  sin  tpi  cos ' 
sin  -ipi  cos  'ipi  sin^  -ipi 


£  s^M 
(Wo/2)  l^i=l 


El  91^ 

c2  •  R2 


COS^  -Ipi 

sin  ipi  cos  'ipi 


sin  ijji  cos  ifji 
sin^  -ipi 


D.2  With  an  Obstacle 

Suppose  that  there  is  an  obstacle  B  at  {Xg,yg)-  This  induces  an  azimuth 
modulation  say  fg{x^,yA-  So  the  attenuation  factors  can  be  written  as, 

Ai{Xj,,yj,)  ^ ^  ^  f bAt^  Vt) 


dAj 

dx^ 


(  Gi\  r  /  dRi{x^,yA  fGi\  dfg{x^,yA 

V  RjjAAT.yT)  \AJ  dx^ 


fg{xT,yT)  cos 'ipi  + 


dfB{x^,yA 

dx^ 
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And  similarly, 


Now, 


^2  J  Jb  sin  Qy^ 


-^j  fs{x^,y^)cos'i/ji- Ri - — - 


f )  fli^r.Vr)  Uosi.,-R,^yX[j^^^ 


—2Ri  cos 'ipi 


d\nfg{x^,y^ 


Similarly,  we  have 


-^i  G-  \ 

^  =  [-A]  fli^T^yr)  sin^^/’i  +  i?; 


Fc^  j  j  B  yT 


-2Ri  sin'^i 


d\nfg{x^,y^) 

dyr 

'dlnfsix^.y^ 

9yT 


dx^J  \dy^ 


^Gl)  f2 
Rf  /  ^ 


fl{x^,y^)  {cosA sin  A+ 


2  ( <91nA(xj,,yj,)\  f  dhifg{x^,y^) 


7?  ,  ( dRifg{x^,y^)\  .  (d\a.fg{xj.,y^) 

-  Ri  cos  il)i  - — 1-  -  Ri  sin  il^i  - f— ' — 1- 

\  ^Ut  )  \  UXrjp 


So,  we  have  the  FIM  as 


I{x^,yT)  -  E*=A  +  ^) /s(A>I/t) 


cos^  -Ipi  sin  Ipi  cos  ipi 


sin  Ipi  cos  Ipi  sin^  ipi 


(D.2) 


where  the  matrix  Jj  is  given  by 


Jjlii  —  Ri 


2  [ dRifp{x^,y^)\  [ d\nfg{x^,y^) 


—  2  Ri  cos  Ipi 
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'i  12  — 


d\nfg{x^,y^)\  fd\nfg{x^,y^) 


dx^ 
-Ri  sin  ipi 


dyj 


d\nfg{x^,y^) 
dx^ 


—  Ri  cos  -ipi 


<91n/g(3;^,i/r) 

QVt 


=  21 


Ji  22  — 


R- 


dlnfsixT^yT) 

dy-r 


2Ri  sin  ipi 


d\nfg{xT,yT) 

dyr 


(D 
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LIST  OF  SYMBOLS,  ABBREVIATIONS,  AND  ACRONYMS 


AENR  Average  energy-to-noise  ratio 
ASNR  Average  signal-to-noise  ratio 
CAM  Complex  cross-ambiguity  matrix 
CRLB  Cramer-Rao  lower  bound 
FDOA  Frequency  difference  of  arrival 
FIM  Fisher  information  matrix 
LPI  Low  probability  of  intercept 
LS  Least  squares 
LSF  Least  squared  error 
MLF  Maximum  likelihood  estimator 
NP  Neyman-Pearson 
PDF  Probability  density  function 
RMS  Root  mean  square 
SNR  Signal-to-noise  ratio 
TDOA  Time  difference  of  arrival 
TOA  Time  of  arrival 
UMP  Uniformly  most  powerful 
WLS  Weighted  least  squares 
WLSE  Weighted  least  squares  error 

{vx,Vy,Vz)  The  velocity  of  the  target  in  the  x,  y  and  2:  directions 
{xi,yi)  Location  of  sensor  i 

{xt,  yTi  zt)  The  x,  y  and  coordinates  of  the  target 
^  Noise  spectral  density 

Tpi  Angle  of  sensor  i  from  the  positive  x-axis  measured  at  the  emitter 
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al  Noise  variance  at  sensor  i 

Ti  Unknown  time  of  arrival  of  signal  at  sensor  i 

A  M  X  1  vector  of  the  unknown  attenuation  factors 

A'  M  —  1  X  1  vector  of  the  unknown  relative  attenuation  factors 

B  Complex  cross-ambiguity  matrix 

Itv  An  N  X  N  identity  matrix 

P  N  X  N  permutation  matrix 

W  N  X  N  DFT  matrix 

A  Complex  attenuation  factor  at  sensor  i 

Ai  Unknown  attenuation  factor  at  sensor  i 

tti,  bi  Fourier  coefficients  of  the  signal 

c  Propogation  speed  of  signal 

dmax  Distance  between  the  farthest  pair  of  sensors 

Fq  Fundamental  frequency  of  the  Fourier  series 

Fg  Sampling  frequency 

ki  Discrete  Doppler  shift  at  sensor  i 

M  Number  of  sensors 

N  Number  of  signal  samples  at  each  sensor. 
rii  Discrete  time  delay  at  sensor  i 
Pd  Probability  of  detection 

PpA  Probability  of  false  alarm 

Ri  Distance  from  the  emitter  to  the  ith  sensor 
Ti  Signal  received  at  sensor  i 

s  Transmitted  Signal 

T  Length  of  observation  interval 

to  Unknown  transmission  time 
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Tg  Non-zero  length  of  the  signal 
Wi  Additive  Gaussian  random  process  at  sensor  i 

Sgi  Signal  energy  at  sensor  i 
X  Information  matrix 

Xgi  Fisher  information  matrix  of  the  unknown  parameter  vector  6 

t]  3x1  vector  of  the  unknown  emitter  location  coordinates  and  the  transmis¬ 
sion  time 

(j)  2N  —  1  X  1  vector  of  the  unknown  Fourier  coefficients 
T  M  X  1  vector  of  the  unknown  TOAs 

t'  M  —  1  X  1  vector  of  the  unknown  TDOAs 

h{t)  2N  —  1  X  1  vector  as  dehned  in  Appendix  B.3 
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